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ABSTRACT 


Fuel  or  battery  consumption  of  unmanned  aerial  vehicles  (UAVs)  can  be 
improved  by  utilizing  or  avoiding  air  currents.  This  thesis  adopts  a  network  modeling 
approach  to  formulate  the  problem  of  finding  minimum  energy  flight  paths.  The  relevant 
airspace  is  divided  into  small  regions  using  a  grid  of  nodes,  inter-connected  by  arcs.  A 
function,  representing  energy  cost,  is  defined  on  every  arc  in  terms  of  the  solution  of  a 
constrained  nonlinear  program  for  the  optimal  local  airspeed  to  fly  in  a  given  wind  field. 
Then,  shortest-path  models  are  implemented  on  the  network  to  find  the  optimal  paths 
from  an  origin  to  a  destination.  Five  models  are  studied  and  they  correspond  to  cases  of 
pre-planning  of  flight  routes  and  dynamic  updating  of  routes  during  the  course  of  the 
flight.  These  models  use  three-dimensional  grids  of  forecasted  wind  currents,  produced 
by  the  Naval  Research  Laboratory’s  Coupled  Ocean-Atmosphere  Mesoscale  Prediction 
System  (COAMPS)  with  horizontal  resolution  of  1  km.  One  of  the  shortest-path  models, 
a  stochastic-dynamic  model,  assumes  real-time  measurement  capabilities  of  the  wind 
velocity  in  the  vicinity  of  the  UAV,  through  its  GPS-INS  system,  and  provides  updated 
waypoints  to  follow  after  every  measurement.  For  each  model,  the  energy  costs  of  the 
shortest-path  solutions  for  1000  randomized  missions  over  a  Nevada  test  site  are 
simulated  and  compared  to  the  energy  costs  of  straight-line  paths.  For  a  100  kg  UAV,  the 
dynamic  model  produces  an  average  reduction  of  15.1%  in  the  energy  consumption  along 
40  km  long  round  trips,  and  an  average  reduction  of  30.1%  under  windy  conditions  with 
average  wind  speeds  larger  than  15  m/s.  A  stochastic-dynamic  model  for  maximum 
duration,  solved  using  a  heuristic  algorithm,  achieves  an  average  increase  of  32.2%  in  the 
flight  duration  for  a  100  kg  UAV. 
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EXECUTIVE  SUMMARY 


The  goal  of  this  study  is  to  provide  a  quantitative  analysis  of  the  potential  energy 
savings  for  unmanned  aerial  vehicles  (UAVs)  through  the  use  of  wind  currents.  The  idea 
is  to  implement  a  shortest-path  network  model,  which  provides  the  minimum-energy 
route  between  an  initial  point  and  a  destination,  using  high-resolution  three-dimensional 
weather  forecasts  and  on-board  dynamic  estimations  of  the  wind  field.  Energy  savings 
can  be  directly  translated  to  fuel  or  battery  savings  and  may  lead  to  longer  flight  distance 
and  duration  as  well  as  improved  operational  flexibility. 

High  resolution  or  mesoscale  weather  forecasts  are  available  from  the  Naval 
Research  Laboratory’s  Coupled  Ocean-Atmosphere  Mesoscale  Prediction  System 
(COAMPS)  model.  This  is  a  state-of-the-art  numerical  model  applicable  to  any  region  of 
the  world  and  provides  72-hours  weather  forecasts  based  on  low  resolution  global 
weather  data  and  the  terrain  features  in  the  region  of  interest. 

On-board  estimations  of  wind  fields  are  available  on  various  UAVs  that  have 
global  positioning  systems  and  inertial  navigation  systems  (GPS-INS  module).  For  such 
UAVs,  wind  velocity  estimates  can  be  achieved  by  subtracting  the  UAV’s  ground 
velocity  vector,  given  by  the  GPS,  from  its  inertial  velocity  vector,  given  by  the  INS. 

The  following  network-based  models  for  finding  minimum-energy  paths  were 
implemented: 

1.  A  deterministic  shortest-path  model  that  assumes  complete  infonnation 
about  the  wind  field  in  the  potential  flight  space.  This  model  provides  an 
upper  bound  for  the  possible  energy  savings. 

2.  A  stochastic-static  shortest-path  model  that  uses  a  forecast  of  wind 
velocities  and  provides  the  path  with  the  minimum  expected  energy  cost. 
This  corresponds  to  the  case  of  a  pre-flight  routing  plan  without  any 
dynamic  updates. 

3.  A  stochastic-dynamic  shortest-path  model  that  assumes  a  given  weather 
forecast,  as  well  as  dynamic  estimations  of  local  wind  fields  through  the 
on-board  GPS-INS  measurements. 


The  flight  space  is  a  three-dimensional  box,  which  is  divided  into  a  grid  of  nodes, 
representing  potential  waypoints,  with  horizontal  spacing  of  1  km  and  vertical  spacing  of 
100  m.  Nodes  are  connected  by  arcs  representing  potential  flight  segments.  Each  node  is 
connected  only  to  its  nearest  neighbors  in  the  general  direction  of  flight  including  those 
along  diagonals. 

Network  representation  of  energy  consumption  is  achieved  through  the 
construction  of  an  aerodynamic  model  of  the  UAV  that  provides  the  arc  costs  for  the 
shortest  path  models.  Since  the  network  is  discretized  into  relatively  short  flight  segments 
represented  by  arcs,  it  is  assumed  that  the  wind  velocity  along  each  arc  is  constant,  and 
that  changes  in  the  wind  field  lead  to  changes  in  the  frame  of  reference  of  the  UAV,  but 
not  to  aerodynamic  instabilities.  Given  a  constant  wind  velocity,  this  study  provides  an 
analytical  fonn  for  the  energy  cost  of  flight  along  an  arc.  This  result  is  given  as  a  function 
of  the  chosen  airspeed,  and  of  the  following  mechanical  and  aerodynamic  parameters  of 
the  UAV: 

•  Mass 

•  Wingspan 

•  Parasite  surface  area 

•  Oswald’s  efficiency  factor 

•  Engine’s  efficiency 

Since  the  UAV’s  airspeed  is  a  free  parameter  that  can  be  changed  through  its 
flight,  this  study  expands  the  arc  cost  formulation  into  a  nonlinear  optimization  model 
that  provides  the  best  speed  to  fly  on  each  arc  to  achieve  minimum  energy  consumption 
on  that  arc.  This  nonlinear  optimization  model  has  the  following  constraints  on  the 
UAV’s  airspeed  and  air-velocity,  for  each  arc  in  the  network: 

•  The  airspeed  must  be  higher  than  the  stall  speed. 

•  The  airspeed  must  be  lower  than  the  maximum  value  that  can  be  achieved 
by  the  engine. 

•  The  air  velocity  must  be  large  enough  to  overcome  the  wind  velocity  and 
achieve  a  flight  on  the  path  defined  by  the  given  arc  with  a  resulting 
ground  speed  that  is  higher  than  a  certain  minimum  value. 


xiv 


The  last  constraint  is  found  through  the  geometry  of  the  problem  by  noticing  that 
the  air-velocity  must  lie  on  the  same  plane  created  by  the  ground  and  wind  velocity 
vectors. 

COAMPS  forecasted  wind  fields  with  horizontal  resolution  of  1  km  are 
transfonned  from  terrain-following  coordinates  (sigma  levels)  into  Cartesian  coordinates 
and  interpolated  on  the  network,  such  that  wind  velocities  are  translated  into  arcs  costs 
through  the  nonlinear  optimization  model.  Vertical  levels  are  linearly  interpolated  into 
100  m  grid  spacing,  in  order  to  provide  constant  small  climbing  and  descending  angles  in 
the  range  of  +5.7  to  -5.7  degrees. 

Following  the  spatial  interpolation  of  the  forecasts,  a  temporal  interpolation  of  a 
few  different  forecasts,  at  different  points  in  time  is  performed.  This  is  done  in  order  to 
compensate  for  the  fact  that  the  UAV  will  reach  different  points  in  space  at  different 
times.  This  pre-processing  interpolation  is  done  under  the  assumption  that  the  UAV  will 
fly  from  origin  to  destination  with  a  representative  speed. 

Wind  speeds  are  modeled  according  to  a  two-parameter  Weibull  probability 
density  function,  and  wind  direction  is  modeled  according  to  a  nonnal  probability  density 
function.  Both  of  these  distribution  functions  are  fitted  in  this  work  using  observed 
quantities  from  the  Marina  weather  station,  California.  This  data  includes  measurements 
of  wind  speed  and  direction  in  altitudes  up  to  1500  m,  throughout  the  year  2002.  The 
vertical  component  of  the  wind  is  usually  insignificant  compared  to  the  horizontal 
components.  Hence,  its  distribution  is  ignored  and,  instead,  the  vertical  forecasted  wind 
speed  is  taken  as  constant. 

Model  I:  Deterministic  Shortest-path 

The  detenninistic  shortest-path  model  uses  statistical  knowledge  about  the  wind 
and  a  wind  forecast  to  create  a  realistic  wind  field  for  the  area  of  interest.  Dijkstra’s 
algorithm  is  then  used  to  find  the  path  with  the  minimum  energy  cost  from  origin  to 
destination.  This  corresponds  to  the  ideal  case,  where  complete  information  about  the  arc 
costs  on  the  network  is  given.  Results  from  this  model  serve  as  an  upper  bound  on  the 
optimal  energy  savings  for  the  other  models. 
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Model  II:  Stochastic-static  Shortest-path 

The  stochastic-static  shortest-path  model  uses  the  statistical  behavior  of  the  wind 
to  calculate  the  expected  (average)  arc  costs  on  the  network  and  use  these  to  find  the 
shortest  path  with  Dijkstra’s  algorithm.  This  model  represents  the  case  of  no  flight  plan 
updates  during  the  flight.  Results  from  this  model  serve  as  a  lower  bound  on  the  energy 
savings  for  the  stochastic-dynamic  model.  The  expected  arc  costs  are  calculated  through 
the  averaging  of  many  realizations,  or  by  the  use  of  a  pre-constructed  table  of  energy 
costs,  for  every  wind  velocity  (one  value  for  every  combination  of  speed  and  direction). 

Model  III:  Stochastic-dynamic  Shortest-path 

The  stochastic-dynamic  shortest-path  model  uses  an  approximate  dynamic 
programming  algorithm  that  provides  the  expected  cost-to-go  from  each  node  to  the 
destination  node.  This  model  dynamically  re-optimizes  the  UAV’s  route  based  on 
dynamically  collected  information  about  the  local  wind  field.  The  calculation  of  the 
expected  cost-to-go  from  each  node  is  carried  out  prior  to  actual  flight  and  only  uses  the 
probability  distributions  of  the  arc  costs.  Specifically,  the  expected  cost-to-go  is 
calculated  through  a  converging  algorithm  that  starts  with  an  initial  guess  for  the 
expected  cost-to-go.  Each  iteration  updates  the  expected  cost-to-go  at  a  node  using  a 
randomly  generated  realization  of  the  on-board  wind  velocity  measurement 
corresponding  to  the  outgoing  arcs  of  the  current  node  as  well  optimization  of  the 
possible  moves  from  this  node.  Given  the  expected  cost-to-go,  the  required  calculations 
during  flight  is  trivial:  At  each  node  in  the  network,  the  UAV  measures  the  wind  and 
determines  the  optimal  arc  to  proceed  on  using  that  measurement  and  the  expected  cost- 
to-go.  This  model  has  a  recursive  nature.  It  assumes  that  the  best  decision  will  be  made  at 
every  point  in  the  future,  based  on  the  infonnation  available  at  that  point. 

Errors  in  the  dynamic  estimation  of  the  wind  velocity  are  simulated  through  the 
use  of  available  data  in  the  literature  regarding  the  accuracy  of  GPS-INS  systems.  A 
“spatial  fluctuation  error”  is  also  taken  into  account  to  represent  the  different  wind 
velocities  in  different  points  along  the  arcs.  This  is  done  by  calculating  the  variation  in 
wind  velocities  for  the  collection  of  all  adjacent  node  pairs  in  the  forecast. 
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Model  IV:  Correlated-based  Stochastic-dynamic  Shortest-path 

This  model  uses  a  correlation-based  stochastic-dynamic  shortest-path  algorithm, 
which  utilizes  the  spatial  correlation  structure  of  the  wind  to  improve  the  accuracy  in  the 
estimation  of  the  distribution  of  wind  velocities.  This  model  assumes  a  spatially- 
autocorrelated  joint  distribution  for  the  wind  velocities,  instead  of  independent  marginal 
distributions  as  was  assumed  in  the  previous  models.  This  model  was  not  implemented  in 
a  simulation,  as  data  on  the  spatial  correlation  structure  of  wind  was  not  available  at  the 
time  of  the  study.  Experimental  data  may  be  achieved  by  simultaneous  measurements  of 
wind  velocities  in  different  altitudes  over  a  set  of  nearby  locations.  The  main  drawback  of 
this  proposed  model  is,  however,  its  expected  heavy  computational  requirements  and 
long  processing  time. 

Model  V:  Stochastic-dynamic  Maximum  Duration  Path 

In  addition  to  the  minimum-energy  models,  a  heuristic  maximum-duration  model 
is  constructed  based  on  the  stochastic-dynamic  model.  This  model  also  assumes  a  back- 
and-forth  flight  path  between  two  points  and  instead  of  minimizing  energy  per  distance,  it 
minimizes  the  energy  per  unit  time,  or  maximizes  the  duration  of  flight  per  unit  of 
energy. 

The  first  three  models  and  the  fifth  model  are  implemented  in  MATLAB,  using 
COAMPS  weather  forecasts  over  a  160  km-by-160  km  region  at  Yucca  Air  Force  Test 
Site,  Nevada.  1000  north-south  and  east- west  paths  of  40  km  round  trips  (approximately 
20  km  ingress  plus  20  km  egress)  are  simulated  through  a  Monte  Carlo  numerical 
simulation.  Origins  and  destinations  points  are  chosen  randomly  within  the  constraint  that 
they  lie  no  less  then  100  meters  above  the  highest  terrain  features  in  the  area  of  flight. 
Each  of  the  models  produces  a  suggested  flight  path  and  an  associated  energy  cost,  which 
is  compared  to  the  energy  cost  of  a  straight  line  flight  path  with  a  constant  speed. 
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As  expected,  the  deterministic  model  provides  an  upper  bound  on  the  amount  of 
energy  savings  for  the  other  models.  The  stochastic-dynamic  model  produces  results  that 
are  significantly  better  than  the  stochastic-static  model.  The  latter  model  produces  energy 
savings  of  no  more  than  10%.  Results  for  the  stochastic-dynamic  model  are  shown  in 
Figures  la  and  lb,  for  representative  10  kg  and  100  kg  UAVs. 

For  a  100  kg  UAV,  the  stochastic-dynamic  model  produces  an  average  reduction 
of  15%  in  the  energy  consumption,  compared  to  the  straight-line  paths,  and  an  average 
reduction  of  30%  under  windy  conditions  with  average  wind  speeds  larger  than  15  m/s. 
The  model  provides  significant  improvements  in  head-wind  scenarios,  in  which  the 
optimized  flight  paths  avoided  stronger  wind  currents.  Conversely,  the  model  does  not 
provide  significant  improvements  in  tail-wind  scenarios,  where  the  optimized  flight  paths 
are  similar  to  straight-line  paths. 
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Figure  1.  Numerical  results  from  the  stochastic-dynamic  model  for  1000  simulations 
of  40  km  round  trips  over  Yucca  Test  Site,  Nevada.  Energy  costs  of  the 
suggested  paths  are  compared  to  the  energy  costs  of  straight-line  paths  for 
(a)  10  kg  UAV  and  (b)  100  kg  UAV. 

The  maximum-duration  model  is  implemented  using  the  same  data  and  1000 
randomly  generated  scenarios,  but  now  without  the  horizontal  degrees  of  freedom.  That 
is,  the  UAV  cannot  deviate  from  a  straight  line,  which  corresponds  to  a  border-patrol 
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scenario.  Results  for  a  100  kg  UAV  are  shown  in  Figure  2.  On  average,  an  increase  of 
32.2%  in  the  flight  duration  is  achieved,  compared  to  a  straight-line  flight  with  a  constant 
speed. 


Figure  2.  Numerical  results  from  the  maximum  duration  stochastic-dynamic  model 
for  1000  simulations  of  40  km  round  trips  over  Yucca  Test  Site,  Nevada. 
Flight  duration  per  unit  energy  was  compared  to  the  same  ratio  in  the  case  of 
straight-line  paths  for  a  100  kg  UAV. 

In  conclusion,  numerical  simulations  has  shown  that  accounting  for  wind  currents 
during  flight  path  planning  may  result  in  a  significant  reduction  of  energy  consumption 
for  various  UAVs  without  any  additional  hardware.  Avoidance  of  strong  head  winds  is 
shown  to  have  the  largest  impact  on  energy  consumption. 


xix 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xx 


ACKNOWLEDGMENTS 


The  author  would  like  to  acknowledge  the  following  contributors  to  the  work: 

Professor  Johannes  O.  Royset,  for  his  thorough  guidance  on  algorithms,  networks 
and  statistical  methods  and  for  constantly  making  sure  that  the  work  is  on  the  right  path,1 

Professor  Kevin  Jones,  for  his  invaluable  advices  on  aerodynamics  and  mechanics 
ofUAVs, 

Dr.  Jason  Nachamkin,  for  his  comprehensive  instruction  on  the  interpretation  and 
analysis  of  COAMPS  weather  forecasts, 

Professor  Richard  (Dick)  Lind,  for  providing  critical  data  for  this  thesis  from  his 
measurements  in  the  Marina  weather  station,  California, 

Maj.  Jim  Wong  and  Ms.  Oh  Pei  Tze  for  their  continuous  help,  encouragement  and 
good  ideas, 

And  the  Naval  Postgraduate  School  and  Naval  Research  Laboratory  staff  who 
shared  their  knowledge  and  insights  with  much  kindness. 


1  Though  not  necessarily  on  the  minimum-energy  path. 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


I.  INTRODUCTION 


A.  MOTIVATION 

Unmanned  aerial  vehicles  (UAVs)  are  expected  to  become  an  important  force 
multiplier  on  the  future  battlefield.  They  will  have  a  profound  impact  in  the  areas  of 
Intelligence,  Surveillance,  and  Reconnaissance  (ISR),  Suppression  of  enemy  defenses 
(SEAD),  Electronic  Warfare  (EW),  and  attack/strike  operations  [1],[2],  The  amount  of 
time  a  UAV  can  stay  in  the  air  and  the  maximum  distance  it  can  reach  are  critical  factors 
for  all  these  missions.  Atmospheric  drag  forces  and  the  UAV’s  energy  resources 
(gasoline  or  battery)  restrict  the  flight  range  of  the  UAVs.  As  UAVs  are  made 
increasingly  smaller,  this  restriction  becomes  more  and  more  prominent  due  to  the 
smaller  weight  allocated  to  battery  and  fuel  compartments.  The  energy  consumption  of  a 
UAV  can  sometimes  be  improved  through  changing  its  structural  design  or  flight  velocity 
[3],  but  also  through  a  better  use  of  available  air  currents  and  gusts.  This  thesis  will 
analyze,  by  means  of  minimum  energy  path  models,  the  possibility  of  using  the  wind  as  a 
way  to  reduce  energy  consumption. 

B.  THESIS  OBJECTIVE  AND  SCOPE 

The  objective  of  this  work  is  to  provide  a  quantitative  analysis  of  the  potential 
energy  savings  for  UAVs  through  the  use  of  wind  currents.  The  idea  is  to  implement 
shortest-path  network  models,  which  provide  the  minimum-energy  or  maximum-flight- 
duration  route  between  an  initial  point  and  a  destination,  using  high-resolution  three- 
dimensional  weather  forecasts  and  on-board  dynamic  estimations  of  the  wind  field.  These 
models  are  based  on  the  stochastic  behavior  of  the  wind  and  utilize  mesoscale  wind 
forecasts  and  on  dynamic  updates  of  the  wind  velocity  by  means  of  onboard 
measurements  throughout  the  UAV’s  flight.  A  significant  portion  of  the  work  focuses  on 
the  conversion  of  aerodynamic  and  physical  characteristics  of  a  flight  in  a  wind  field  into 
a  network  model  formulation. 
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The  proposed  models  are  capable  of  handling  complex  wind  structures.  At  low 
altitudes  -  up  to  2  km  -  there  are  usually  more  frequent  changes  in  the  wind  direction. 
Conversely,  at  higher  altitudes,2  wind  currents  are  relatively  constant  and  the  optimal 
flight  paths  are  usually  easy  to  predict.  Hence,  this  thesis  targets  UAV  missions  at  low 
altitudes  where  route  optimization  models  may  be  most  beneficial.  Since  weather 
forecasts  are  currently  available  at  the  horizontal  resolution  of  1  km,  this  study  focuses  on 
mission  with  flight  lengths  significantly  greater  than  1  km. 

C.  POTENTIAL  BENEFITS  OF  THE  STUDY 

Energy  conservation  can  be  directly  translated  to  fuel  and  battery  savings  and  may 
lead  to  longer  flight  distance  and  duration,  and  a  higher  flexibility  in  real-time  decisions 
regarding  a  UAV  mission.  Implementation  of  the  proposed  shortest-path  models  requires 
only  software  modifications  in  the  control  station  of  the  UAV.  No  hardware  changes  are 
typically  needed  as  most  UAVs  already  have  the  necessary  instruments  such  as  global 
positioning  system  (GPS)  and  inertial  navigation  system  (INS)  for  measuring  wind. 

D.  RELATED  WORK 

There  are  several  studies  on  optimal  routing  of  UAVs  with  respect  to  energy 
resources,  target  acquisition  and  obstacles  avoidance.  These  studies  have  considered  both 
deterministic  and  stochastic  situations  and  have  derived  methods  based  on,  for  example, 
potential  fields  [4],  graph  search  methods  [5]  and  evolutionary  algorithms  [6], [7], [8]. 

Paul  McCready  published  his  theory  in  the  early  1950s  on  rules  for  efficient 
gliding,  and  specifically,  guidelines  for  optimal  flight  between  thennals  (“MacCready 
Speed,”  see  [9]).  Since  then,  these  guidelines  have  served  as  a  significant  help  to  pilots’ 
intuition.  A  large  amount  of  work  has  since  been  done  on  energy  efficient  speeds  and 
trajectories  for  gliders  and  utilization  of  thennals  for  “soaring.”  These  studies  include  the 
use  of  optimal  control,  parameter  optimization  and  the  variational  fonnulation  [9] -[13],  to 
achieve  improved  paths  (higher  or  faster  flight)  using  local  observations  of  vertical 
currents.  These  models  can  be  helpful  in  UAV  missions  either  through  reduction  of 


2  See  Figure  11,  p.  17. 
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energy  consumption  by  leveraging  the  wind,  or  by  reducing  the  engine’s  noise  for  the 
purpose  of  concealment  [14].  Due  to  the  unpredictable  behavior  of  thermals,  the  models 
in  these  studies  provide  only  locally-optimal  paths  under  initial  constraints,  and  do  not 
make  use  of  large  scale  wind  currents  forecasts. 

Large-scale  optimization  of  flight  duration  and  fuel  consumption  is  frequently 
used  in  commercial  aircrafts.  This  is  done  through  the  use  of  jet  streams  and  real-time 
weather  forecasts  of  high-altitude  wind  currents  [15].  The  stochastic  and  dynamic 
behavior  of  weather  obstacles,  such  as  storms  is  also  studied  [16].  In  the  world  of  UAVs, 
evolutionary  algorithms  are  used  to  find  an  adaptive  optimized  path,  using  large-scale 
wind  fields,  e.g.,  for  a  Trans-Pacific  Crossing  [8],  or  adaptive  obstacle  avoidance  [6], 

None  of  the  above  studies  utilize  low-altitude,  high-resolution  weather  forecasts 
for  minimum-energy  or  maximum-duration  paths.  This  thesis  will  use  such  forecasts  and 
formulate  network  models  that  provide  optimal  paths  for  the  case  of  low-altitude  UAV 
flights. 

E.  THESIS  ORGANIZATION 

Chapter  II  discusses  the  nature  of  wind,  presents  an  analysis  of  weather  data  from 
the  Marina  weather  station,  and  describes  COAMPS  mesoscale  weather  forecast  model  of 
the  Naval  Research  Laboratory  (NRL).  Chapter  III  discusses  the  construction  of  a 
network  model  and  an  aerodynamic  model  for  the  energy  cost  of  flight  along  an  arc  in  the 
network.  Chapter  IV  describes  the  five  network  models  for  minimum-energy  paths:  a 
deterministic  model  that  assumes  perfect  knowledge,  a  stochastic-static  model  for  a  pre¬ 
flight  plan,  a  stochastic-dynamic  model  that  takes  into  account  dynamic  measurements,  a 
correlation-based  stochastic-dynamic  model  that  considers  the  correlation  structure  of  the 
atmosphere  and  a  maximum-duration  model  for  the  maximization  of  the  flight  time. 
Chapter  IV  discusses  the  implementation  of  the  models  in  a  simulation  of  flights  above  a 
Nevada  test  site,  and  presents  numerical  results.  Chapter  VI  concludes  the  work. 
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II.  WIND  AND  FORECASTS 


A.  CHAPTER  OVERVIEW 

Our  study  is  motivated  by  the  potential  of  energy  savings  for  unmanned  aerial 
vehicles  (UAVs)  through  the  use  of  wind  currents.  In  this  chapter,  the  statistical 
characteristics  of  wind  are  discussed,  followed  by  analysis  of  weather  data  collected  in 
the  Marina  weather  station  and  a  description  of  COAMPS  weather  forecasts.  The  general 
way  in  which  wind  data  and  forecasts  will  be  implemented  in  our  network  models  is  also 
discussed.  This  includes  the  temporal  and  spatial  correlation  structure  of  the  wind 
velocity,  which  serves  as  a  basis  for  a  correlation-based  shortest  path  model. 

B.  WIND  BEHAVIOR 

Wind  -  the  flow  of  air  in  the  atmosphere  -  is  caused  by  pressure  gradients.  On  the 
global  scale,  the  major  driving  factors  for  wind  are  the  sun’s  uneven  heating  of  the  earth 
and  the  Coriolis  force,  which  is  created  by  the  rotation  of  the  earth.  In  mesoscale  wind, 
which  is  on  the  scale  of  hundreds  of  meters  up  to  tens  of  kilometers,  the  main  driving 
effects  are  differential  heating3  and  terrain  features.  Micro-scale  winds  exist  on  the  scale 
of  tens  to  hundreds  of  meters  (e.g.,  thennals),  and  are  essentially  unpredictable. 

On  small  scales,  wind  speed  and  wind  direction  can  both  be  seen  as  random 
variables.  The  observed  quantities  will  be  distributed  around  some  mean  value.  Wind 
speeds  are  frequently  modeled  according  to  a  two-parameter  Weibull  probability  density 
function  [  1 7]-[20] .  This  thesis  assumes  that  the  wind  direction  has  a  nonnal  probability 
density  function  with  the  forecasted  direction  as  its  mean.  Both  of  these  distribution 
functions  are  fitted  in  this  work  using  observed  quantities. 


3  For  example,  in  the  case  of  a  sea  breeze,  the  source  for  differential  heating  is  the  land’s  faster  heating 
during  the  day  and  cooling  during  the  night  compared  to  the  sea,  which  leads  to  a  lower  pressure  above  the 
sea  during  the  days  and  a  higher  pressure  throughout  the  evenings  and  nights. 
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The  Weibull  distribution  has  the  following  probability  density  function: 


(1)  f(V)  = 


V  >  0,  k>0,  c>0 . 


and  the  cumulative  probability  distribution  function: 


(2)  F(V)=r  f(V)dV=  1-exp  -  - 

J  —co  1 


where  V  is  the  wind  speed  and  k  and  c  are  the  two  parameters  of  the  distribution.  The 
mean  and  variance  of  the  Weibull  distribution  are  given  by: 


(3)  E{V)  =  cT\\  +  - 


Var(V)  =  c2  r  1  + - r  1  +  — 


where  T  is  the  gamma  function  in  the  positive  domain: 

oo 

(5)  T(z)  =  jV”1  exp(-t) dt,  z>0. 


There  are  several  methods  of  determining  the  parameters  c  and  k  .  One  of  them 
[18]  uses  the  mean  and  standard  deviation  of  the  wind  speed  to  find  c  and  k  as  follows: 


(6)  k=  ^ 

{  V 


r  i+ 
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Here,  V  is  the  mean  wind  speed,  and  ov  is  the  standard  deviation  of  the  wind  speed. 
This  method  is  used  here  together  with  wind  observations  and  the  forecasted  wind  speed 
values,  as  detailed  in  section  II. C. 

A  random  variable  V  that  follows  a  Weibull  distribution  can  be  created  using  a 
uniformly  distributed  random  variable  and  the  Inverse  Method  for  the  cumulative 
distribution  function  in  eq.  (2).  In  this  work,  uniformly  distributed  random 
variables  /  ~  Un{ 0,1)  are  generated  using  MATLAB’s  random  variable  generator  and  the 
random  wind  speed,  V  ,  is  calculated  through: 

(8)  V  =  F  l{%)  =  [—  ln(l  -  x)]Uk  c . 


The  normal  distribution  of  the  wind  direction,  6 ,  has  the  following  parameters: 


(9)  g(0)  =  - - exp 


'2.71(7 a 


r  [o-o«r 

2<J  n2 


where  og  is  the  standard  deviation  and  #0  is  the  average  wind  direction.  The  random 

variable#  is  created  using  MATLAB’s  generator  for  normally  distributed  random 
variables.  It  is  assumed  that  the  forecasted  wind  direction  is  an  unbiased  prediction  and 
equals  its  average  value  #0 .  <jg  is  calculated  using  observation  data  and  forecasted  wind 
speed  values  as  detailed  in  the  following  section. 


C.  WIND  DATA 


In  the  proposed  models,  the  parameters  of  the  two  distribution  functions  for  the 
speed  and  direction  are  computed  from  data  from  the  Marina  weather  station,  California, 
for  the  year  2002,  collected  by  the  Department  of  Meteorology  of  the  Naval  Postgraduate 
School  (NPS)  [21].  This  data  was  collected  continuously  at  40  different  altitudes  up  to 
1500  m  and  averaged  over  periods  of  30  minutes  using  NPS  915  MHZ  wind  profiler,  a 
Doppler  based  three  beams  radar  for  measurements  of  the  three  components  of  the  wind 
[22], [23], 
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Analysis  of  the  data  collected  reveals  a  consistent  relation  between  crv  and  V ,  and 

also  a  consistent  relation  between  cre  and  V  .  These  relations  are  used  as  a  basis  for  our 
simulations.  The  following  approximated  relations  are  found  through  fitting  of  the 
observed  relations  by  a  polynomial  of  the  third  degree  for  V  : 

ay  f-5.6-10^F+1.75-10~2F2-0.178F  +  0.701 ,  0 < F <  15 

OW  77  ~  }  —  ’ 

v  [0.1,  V>\5 

/11X  f-5  -10  2F3  +  1.63F2 -17.63K  + 69.41  ,  0  <  F  <  15 

(H)  °g*\  - 

6,  V>\5 


Here,  ay  and  V  are  given  in  m/s,  and  ay  is  given  in  degrees.  These  approximations  for 
the  statistical  measures  of  the  wind  are  found  by  dividing  the  year  of  observations  into 
two-hour  periods  in  order  to  achieve  results  for  temporally  local  fluctuations  while 
allowing  enough  observations  (of  30  minutes  intervals)  per  calculation. 


Comparison  of  these  approximations  to  the  observed  data  is  shown  in  Figures  3 
and  4.  Less  data  is  available  for  higher  wind  speeds  (>15  m/s).  Hence,  a  constant  fitted 
value  is  used  to  avoid  over- fitting  of  the  fluctuations. 


Quantification  of  Stochastic  Speed  Measures 
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Figure  3.  Plot  of  av/V  as  a  function  of  V  for  data  collected  at  Marina  Weather 
Station  in  2002  with  altitudes  up  to  1500  m. 
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Quantification  of  Stochastic  Direction  Measures 
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Figure  4.  Plot  of  ag  as  a  function  of  V  for  data  collected  at  Marina  Weather  Station 
in  2002  with  altitudes  up  to  1500  m. 

Eqs.  (10)  and  (11)  are  implemented  in  this  study  using  V  as  the  forecasted  wind 
speed,  and  0O  as  the  forecasted  wind  direction,  such  that  the  standard  deviations  <7V  and 

<jg  are  calculated  from  V  only. 


D.  MESOSCALE  WEATHER  FORECASTS 

Mesoscale  weather  forecasts  are  characterized  by  horizontal  grid  resolutions  on 
the  scale  of  ten  kilometer  down  to  hundreds  of  meters.  These  forecasts  are  created  by 
modeling  the  atmosphere  of  a  certain  region  as  a  three-dimensional  array  of  cells,  where 
each  cell  is  characterized  by  its  physical  attributes,  such  as  air  temperature,  pressure, 
density,  humidity  etc.  In  these  forecast  models,  application  of  numerical  models  for  the 
governing  physical  equations  of  the  atmosphere,  given  initial  and  boundary  conditions, 
leads  to  a  set  of  forecasts  for  a  collection  of  points  in  time. 

COAMPS-Coupled  Ocean-Atmosphere  Mesoscale  Prediction  System-is  a  state- 
of-the-art  weather  forecast  model  for  high  resolution  short-term  (up  to  72  hours) 
predictions  built  by  Naval  Research  Laboratory  (NRL)  [26].  This  model  provides  various 
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weather  characteristics  in  time  intervals  of  down  to  one  hour  and  with  horizontal  spatial 
resolution  of  down  to  1  km.  Vertical  levels  are  given  in  the  fonn  of  sigma  levels,  or 
terrain  following  coordinate  system  [26].  In  order  to  convert  the  forecasts  into  Cartesian 
coordinates,  the  definition  of  the  terrain-following  coordinates  is  used: 

(12)  (7  =  s^^-, 

s-zG 


or: 


(13) 


z  =  zG+  <j 


f  z  3 


V 


s  J 


where  cr  is  the  sigma  coordinate,  zG  is  the  ground  altitude,  5  is  the  maximal  altitude  in 

the  model  and  z  is  the  desired  Cartesian  altitude.  COAMPS  uses  60  sigma  levels,  for 
altitudes  up  to  30  km,  as  shown  in  Figure  5. 


Figure  5.  Decreasing  vertical  sigma  levels,  with  respect  to  the  ground  altitude,  in 
terrain-following  coordinates,  from  top  to  bottom. 

Various  tests  were  conducted  to  assess  the  accuracy  of  COAMPS  forecasts,  see 
[23], [24], [27].  Measured  values  of  wind  speeds  were  recently  compared  to  the  forecasted 
values  by  Nachamkin  et  ah,  [27].  Measurements  were  made  in  10  m  altitude,  and 
produced  the  results  of  Table  1  for  low  wind  speeds  (up  to  4  m/s),  with  a  1  km  grid 
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spacing  model.  These  measurements  show  that  the  biases  in  the  forecasted  wind 
measures  are  significantly  smaller  than  the  standard  deviations  and  hence,  this  thesis 
assumes  that  the  wind  forecasts  are  unbiased  estimators  for  the  distributions  of  wind 
speeds  and  directions. 

Table  1 .  Validation  of  COAMPS  predicted  wind  speeds  (From:  [27]). 


Speed  bias 

0.53  m/s 

Speed  standard  deviation 

1.75  m/s 

Direction  bias 

1.28  degrees 

Direction  standard  deviation 

75.2  degrees 

E.  USING  FORECASTED  WIND  FIELDS 

This  study  uses  only  the  three  components  of  the  wind  velocity  in  every  point  in 
the  grid,  and  does  not  make  use  of  the  various  other  physical  parameters  provided  by  the 
forecasts  [26].  Let  us  denote  these  components  by:  Vu(x,y,z,t) ,  Vv(x,y,z,t )  and 

Vw(x,y,z,t ),  where  x ,  v  and  z  are  the  spatial  coordinates  in  the  region  of  interest,  t  is 

the  point  in  time  for  which  the  forecast  refers  to  and  u ,  v  and  w  denote  each  of  the 
Cartesian  wind  components  in  the  x ,  y  and  z  directions,  respectively. 

After  performing  a  spatial  transformation  of  the  terrain-following  coordinates  into 
a  Cartesian  grid,  through  eq.  (13),  wind  components  are  interpolated  into  a  constant 
vertical  spacing,  using  a  linear  approximation.  This  is  done  in  order  to  assure  low  angles 
of  climb  and  descent  when  implementing  a  simulated  flight  and  avoid  non-linearity  in  the 
drag-to-lift  ratio  [29].  Since  the  differences  between  forecast  times  are  small  (down  to 
one  hour),  linear  approximation  for  the  change  in  the  three  components  of  the  wind 

velocity  is  used  in  the  interpolations.  Differences  in  the  forecasted  wind  speed,  V wind 

and  wind  direction,  0  ,  within  3  hours,  are  presented  as  contour  plots  over  a  50  km  by  50 
km  area  in  Figures  6  and  7.  Angles  are  taken  with  respect  to  the  north,  counterclockwise. 
It  is  noticed  that  the  wind  field  experiences  a  general  acceleration  and  a  clockwise  shift  in 
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the  later  forecast,  compared  to  the  earlier  one.  The  three  components  of  the  wind  are 
shown  for  different  altitudes  in  Figure  8.  Generally,  the  vertical  component  is  negligible 


compared  to  the  horizontal  ones. 


9  [degrees] 


Figure  6.  Wind  speed  and  direction  fields  for  two  forecasts,  at  6  am  and  9  am  over  a 
Yucca,  Nevada  test  site,  at  altitude  1000  m. 


Difference  in  Wind  directions  Difference  in  Wind  Speeds 


(a) 


(b) 


Figure  7.  Histograms  of  (a)  difference  in  wind  directions  and  (b)  difference  in  wind 
speeds,  for  three-hours-separated  COAMPS  forecasts  at  6  am,  and  9  am 
over  Yucca,  Nevada  test  site,  at  altitude  1000  m. 
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vd  [m/s] 

50  km  x  50  km  Forecasted  Wind  Field 


582  m  1596  m  2814  m 


Figure  8.  Three  components  of  the  wind  for  weather  forecasts  produced  by  COAMPS 
for  three  different  altitudes  (left  to  right):  582  m,  1596  m  and  2814  m,  over 
Yucca,  Nevada  test  site. 


Since  this  study  uses  a  stationary  probability  model  of  the  wind  it  will  not  matter 
when  the  UAV  arrives  at  an  arc-the  arc  will  always  have  the  same  wind  velocity 
distribution.  Therefore,  a  temporal  interpolation  of  a  few  different  forecasts,  at  different 
points  in  time  follows  the  spatial  interpolation  of  the  forecasts.  This  is  done  in  order  to 
compensate  for  the  fact  that  the  UAV  reaches  different  points  in  space  at  different  times. 
This  pre-processing  interpolation  is  done  under  the  assumption  that  the  UAV  will  fly 
from  origin  to  destination  with  a  representative  speed. 
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Temporal  interpolation  of  forecasts  is  done  through  the  following  equation: 

(14)  V*d(r )  =  V*d  (x,y,z)  =  Vld(x,y,z,tl)  +  a[  V2  d  (x,  y,  z,t2)~  Vxd  (x,  y,  z,tx)). 


where  r  =  (x,y,z)is  a  coordinate  of  V* d  ,  the  interpolated  grid  component  d  e  [u,  v,  vv} , 

Vld  is  the  forecast  at  time  t] ,  V2d  is  the  forecast  at  time  t2  and  a  is  the  interpolation 
factor,  calculated  by: 


(15)  a(x,y,z)  = 


r  -  rn 


t2  t\ 


_  V(^-^o)2  +{y-y0)2  +{z-zo)2 

Vo  (^2  —  h) 


where  ro  =  (x0,  v0,z0)  is  the  coordinate  of  the  origin  and  V0  is  an  estimation  for  the 
UAV’s  speed  along  its  flight  path. 

It  is  noticed  that  for  an  interpolation  at  a  point  in  time  between  two  forecasts  to  be 
possible,  and  not  become  an  extrapolation,  a(x,y,z )  must  obey  the  following  equation: 

(16)  0<a(x,y,z)<l. 


Projection  of  the  grid  on  the  Earth’s  coordinates  is  unnecessary  in  our  models, 
since  the  scales  of  interest  are  of  dozens  of  kilometers,  in  which  the  curvature  of  the  earth 
is  not  significant. 


F.  TEMPORAL  AND  SPATIAL  CORRELATION  OF  THE  ATMOSPHERE 


Temporal  and  spatial  correlation  characteristics  of  the  atmosphere  may  be 
incorporated  in  a  minimum-energy  model  for  better  predictability  of  the  wind  behavior 
and  improved  energy  conservation.  This  section  discusses  the  correlation  structure  of  the 
wind,  in  the  context  of  shortest-path  models. 

Since  forecasted  wind  fields  are  assumed  to  have  the  average  values  of  the  three 
components  of  the  wind  distribution,  this  study  is  interested  in  the  structural  correlation 
of  the  fluctuation  of  the  wind  around  its  average  values.  Correlation  may  be  calculated 
for  any  two  points  in  time  and  space.  Models  for  the  complete  correlation  structure  of  the 
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atmosphere,  in  general,  or  experimental  data  for  specific  regions  was  not  available  at  the 
time  of  the  study.  However,  analysis  of  the  spatial  correlation  at  different  altitudes  and 
the  temporal  correlation  at  one  location  reveals  partially,  as  seen  in  Figures  9  and  10,  the 
complete  space  and  time  correlation  structure  of  the  wind  fluctuations. 


Correlation  as  a  Function  of  Vertical  Distance 
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Figure  9.  Spatial  correlation  of  the  fluctuations  in  wind  speed  at  increasing  altitudes 
with  respect  to  the  wind  speed  at  the  lowest  altitude  available  in  five 
weather  stations  in  California. 


It  is  observed  in  Figure  9  that  the  correlation  in  the  wind  speed  remains  above  0.8 
for  spatial  separations  of  up  to  300  m.  Figure  10  reveals  that  also  the  temporal  correlation 
of  the  wind  speed  and  direction  remains  quite  high  for  the  short  term.  In  these 
experiments,  a  correlation  higher  than  0.8  is  observed  for  up  to  70  minutes  regarding  the 
wind  speed,  and  up  to  25  minutes  regarding  the  wind  direction.  High  correlation  values 
may  be  used  for  better  estimation  of  the  wind  at  different  locations  and  times  given  local 
measurements,  and  therefore  provide  better  results  for  the  shortest  path  models. 
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Figure  10  presents  measurements  of  the  wind  speeds  and  directions  at  two-minute 
intervals.  The  temporal  correlation  structure  is  calculated  according  to  these  periodic 
measurements  recorded  throughout  2002.  High  correlation  is  seen  at  intervals  of  exactly 
24  hours  due  to  the  consistent  weather  cycle  near  the  ocean  (e.g.,  sea  breezes  in  the 
evenings.) 


Temporal  Pattern  of  the  Horizontal  Wind  Direction 


Autocorrelation  of  Horizontal  Wind  Direction  as  a  Function  of  Time  Difference 


(c) 


(d) 


Figure  10.  (a)  Example  of  the  wind  speed  pattern  throughout  three  days  of 

measurements  in  Marina  Weather  station  10  m  above  ground  level,  (b) 
Temporal  correlation  of  the  wind  speed,  (c)  Horizontal  wind  direction 
pattern  for  the  same  case,  (d)  Temporal  correlation  of  the  wind  direction. 
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G. 


WIND  FIELD  COMPLEXITY 


The  shortest-path  models  were  tested  for  flights  at  different  altitudes.  It  was 
observed  that  when  the  wind  field  is  approximately  constant,  then  the  resulting 
minimum-energy  paths  became  similar  to  the  straight  line  paths.  When  the  simulated 
flight  altitude  increases,  a  similar  result  is  noticed.  Therefore,  the  effect  of  altitude  on  the 
complexity  of  the  wind  field,  or  the  amount  of  variation  in  the  wind  speed  and  direction 
as  a  function  of  the  altitude,  is  analyzed  in  Figure  11.  Generally,  the  standard  deviation  of 
the  wind  directions  at  the  same  altitude  is  low  for  altitudes  of  2  km  and  above.  Variation 
in  the  speed  reaches  its  highest  value  just  below  2  km.  Thus,  it  is  assumed  that  the  models 
presented  in  this  work  will  provide  a  significant  improvement  mostly  for  low-altitude 
flights  below  2  km. 


Variations  in  the  Wind  Field  for  different  Altitudes 


(a) 


(b) 


Figure  11.  Observed  variations  in  the  wind  field  as  a  function  of  altitude,  up  to  10  km, 
according  to  a  160  km  x  160  km  COAMPS  forecast  over  Yucca,  Nevada 
test  site,  for  (a)  wind  speed  and  (b)  wind  direction. 
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III.  NETWORK  MODELING  OF  THE  PROBLEM 


A.  CHAPTER  OVERVIEW 

In  this  chapter,  a  network  model  of  the  airspace  is  constructed,  as  well  as  an 
aerodynamic  model  for  the  energy  cost  of  a  flight  on  each  arc  in  the  network. 

B.  NETWORK  MODEL 

The  problem  of  finding  a  minimum-energy  path  between  an  initial  UAV  location 
and  a  destination  can  be  solved  using  network-based  shortest-path  models  and 
algorithms.  The  airspace  is  divided  into  a  grid  of  nodes,  which  represent  potential 
waypoints.  Adjacent  nodes  are  connected  by  arcs,  which  represent  potential  flight  paths 
between  neighboring  waypoints.  Let  us  define  a  network  (N,E) ,  where  N  is  the  set  of 
nodes  and  E  is  the  set  of  arcs.  Arcs  are  denoted  by  ( i,j)eE ,  where  i,jeN  are 
adjacent  nodes. 

C.  NETWORK  ARCHITECTURE 

The  spatial  orientation  of  the  network  is  defined  such  that  the  grid  is  aligned  with 
the  direction  of  the  straight-line  between  the  initial  UAV  location  and  its  destination.  The 
network  can  then  be  built  around  the  straight-line  path,  and  the  minimum-energy  path  can 
be  directly  compared  to  the  straight-line  path. 

In  order  to  contain  the  network’s  size  and  subsequently  the  runtime  of  the 
algorithms,  the  flight  is  restricted  to  always  have  a  “forward”  component.  In  other  words, 
the  only  outgoing  arcs  from  each  node  will  be  those  that  are  directed  towards  the 
neighboring  nodes  in  the  forward  plane,  as  depicted  in  Figure  12.  Every  node  that  is  not 
on  the  boundary  of  the  network  has  nine  outgoing  arcs  in  the  general  direction  of  flight. 
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General  Flight  Direction 


Figure  12.  Example  of  outgoing  arcs  from  five  nodes  in  the  network. 

The  horizontal  spacing  is  1  km,  whereas  the  vertical  spacing  is  100  m.  This  leads 
to  four  potential  arc  lengths:  1  km,  1.005  km,  1.41  km,  1.42  km  and  produces  small 
angles  of  climb  and  descent.  As  mentioned  in  Chapter  II,  these  small  angles  help  in 
avoiding  non-linearity  in  the  drag-to-lift  ratio  [29]  and  maintains  the  validity  of  the 
aerodynamic  expressions. 

Nodes  and  arcs  are  represented  in  a  forward-reverse  star  data  structures  to 
facilitate  quick  access  of  information  by  the  different  algorithms  [30].  In  order  to  avoid 
individual  consideration  of  every  type  of  nodes  (i.e.,  comers/edges/inner  nodes  -  27  types 
in  total  for  a  three  dimensional  rectangular  network)  during  the  initialization  of  the 
networks,  none  of  the  boundary  nodes  contain  any  outgoing  arcs  to  the  rest  of  the 
network.  This  ensures  that  the  inner  nodes  will  still  have  exactly  nine  outgoing  arcs  to  the 
forward  layer,  while  reducing  extra  checks  for  those  nodes  that  have  between  one  and 
eight  outgoing  arcs  to  the  forward  layer.  If  this  individual  “care”  was  given  for  every 
node  type  in  the  network  throughout  the  execution  of  the  shortest-path  algorithms  (e.g., 
Dijkstra’s  [30]),  it  would  result,  for  example,  in  a  computational  cost  of  O(n')  for  a 
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cubic  network,  where  n  is  the  number  of  nodes  in  one  axis  of  the  cube.  The  price  is  an 
extra  0{n 2)  for  every  execution  of  the  shortest-path  algorithm,  since  all  boundary  nodes 
are  practically  meaningless  to  the  solution,  but  are  still  taken  into  account  since  arcs  are 
directed  towards  them,  in  general.  This  turns  out  to  be  useful  particularly  in  cases  where 
statistics  are  made  using  many  different  path  samples  of  a  certain  area,  and  the  algorithm 
is  executed  once  for  every  initialization  of  the  network,  so  the  extra  0(n  ' )  operations  are 
essentially  replaced  by  extra  0(n2)  operations. 


D.  ARC  OPTIMIZATION 


Before  discussing  the  network-scale  shortest-path  optimization,  this  study 
investigates  the  optimal  way  to  cross  an  arbitrary  arc  in  the  network,  under  the  UAV’s 
two  degrees  of  freedom:  flight  speed  and  direction.  This  optimization  is  referred  to  as 
“arc  optimization.” 


1.  Flight  Speed  Optimization  in  Zero  Wind 


In  [3],  Carson  shows  that  the  most  fuel-efficient  flight  speed  (later  known  as 
“Carson’s  speed”),  in  the  case  of  zero  wind,  can  be  expressed  as  a  function  of  the 
aircraft’s  aerodynamic  characteristics.  This  study  will  follow  Carson’s  framework  [3], 
and  assume  a  drag-to-lift  ratio  that  can  be  approximated  as: 

(17)  —  =  AV2  +B/V2, 

L 


where  V  is  the  flight  speed  and,  A  and  B  are  given  by: 
P(h)f 


(18)  A  = 


2  W 


2  W 

(19)  B  = 

p(h)b  ne 


Here,  W  =  mg  is  the  aircraft’s  weight,  given  by  multiplication  of  its  mass  by  the 

gravitational  acceleration,  p(h)  is  the  air  density  at  altitude  h  ,  b  is  the  wing  span,  /  is 

the  parasite  area  of  the  aircraft  and  e  is  Oswald’s  efficiency  factor  of  the  aircraft. 
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This  formulation  deviates  from  the  standard  convention,  which  represents  the 
aerodynamic  parameters  in  a  coefficient  form.  This  representation  was  adopted  in  order 
to  allow  us  to  analyze  the  dependence  of  the  results  on  the  various  aircraft’s  parameters. 
Carson’s  work  was  published  before  the  UAV  era,  and  refers  to  small  manned  aircraft 
only.  An  assumption  in  this  study  is  that  the  above  lift  to  drag  ratio  will  still  hold  for  most 
fixed-wing  UAV’s,  due  to  their  similar  general  structure  to  small  manned  aircraft. 

Equations  (17)  through  (19)  may  not  be  good  approximations  for  the  real 
dependence  of  drag-to-lift  on  speed  when  the  Reynolds  number  is  small,  or  when  the 
UAV’s  structure  is  irregular.  This  problem  is  avoided  in  the  case  of  constant  airspeed, 
where  D  /  L  is  constant  and  results  are  presented  in  the  form  of  the  relative  amount  of 
fuel/battery  capacity  that  is  saved.  This  is  a  result  of  the  fact  that  energy  consumption  can 
be  shown  to  be  proportional  to  the  D  /  L  ratio,  with  or  without  the  presence  of  wind. 

In  cases  where  the  speed-dependent  drag-to-lift  ratio  of  the  UAVs  can  be 
expressed  in  a  more  accurate  form  (e.g.,  based  on  experiments),  this  form  should  replace 
the  generalized  ratio  in  equation  (17)  to  achieve  more  realistic  results. 

The  optimal  lift-to-drag  ratio  occurs  when  eq.  (17)  reaches  its  minimum  value.  If 
the  dependence  of  the  engine’s  efficiency  on  the  aircraft’s  speed  is  ignored,  the  flight 
speed  that  produces  the  optimal  lift  to  drag  ratio  becomes  [3] 


(20)  Vopt  = 


B_ 


=  (2  W)X/2(p2b2fne)l/4, 


and  the  optimal  lift  to  drag  ratio  is: 


(21) 


rD 


\opt 

J 


2y[AB  =2, 


/ 

nb2e 


In  the  arc  optimization  analysis,  where  the  optimal  flight  within  individual  arcs  is 
found,  the  objective  function  is  defined  as  the  total  amount  of  energy  spent  over  the  given 
arc.  This  objective  function  is  measured  in  Joules  as  well  as  the  objective  function  of  the 
network-scale  shortest  path  optimization.  Both  of  them  need  to  be  minimized,  starting 
with  the  arc  optimization. 
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For  a  zero-wind  scenario,  the  total  energy  spent,  or  the  cost  of  crossing  an  arc 
(/',/)  assuming  constant  speed  and  steady  wind  along  the  arc  is  given  by  eq.  (22). 


(22) 


cost{j  = 


Energy  spent  to  resist  air  drag  Air  Drag  Force  ■  Distance  flown 


Engine's  efficiency  Engine's  efficiency 

where  the  engine’s  efficiency  is  defined  here  for  a  constant  speed  per  unit  time: 

c  .  _  _  Energy  dissipated  to  the  air  through  direct  drag 

v  )  hngine  s  ejjiciency  =  ije ngine  —  ~  :  :  7  :  :  ~ 

Potential  energy  in  jaet  or  battery 


The  last  equality  in  eq.  (22)  emerges  from  the  assumption  of  constant  wind,  i.e.,  constant 
drag  force  D  along  the  arc  (/,/) ,  with  length  Xtj : 

(24)  J  D(x)dx  =  DXy . 

arc(i,j) 


All  of  the  processes  that  lead  to  loss  of  energy,  including  drag,  will  eventually  be 
translated  to  heating  of  the  environment.  Air  drag  is  the  main  process  of  energy  transfer 
from  the  engine  to  the  outside  world  for  reasonably  efficient  UAVs.  Hence,  this  study 
introduces  the  engine’s  efficiency  as  a  parameter  that  compensates  for  the  non  drag- 
related  processes. 

Both  the  drag  forces  and  the  engine’s  efficiency  are  functions  of  the  airspeed. 
Using  the  fact  that  the  aircraft’s  weight  equals  its  lift,  the  energy  cost  can  be  expressed 
as: 


(25)  cost ‘  =  DXijrjt 


_1.  =  WX 

engine  ij 


'  D' 

,~L. 


V'enzine  =  ™gX ii 


'  D^ 

J, 


He 


Assuming  constant  engine  efficiency,  for  simplicity,  and  using  Carson’s  optimal 
speed  from  equation  (20)  yields: 


(26)  cost  if1  =  WXU 


f  D 


K  L  hv0„) 


f 

71  1  =  2  WX  ^  77  X 

I  engine  v  ij  ^  (f  en&ne  * 
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2. 


Effect  of  Wind 


Up  to  this  point,  Carson’s  formulation  was  used  for  the  case  of  zero  wind.  This 
formulation  can  be  easily  transformed  to  the  case  of  one-dimensional  wind,  by  adding  the 
wind  velocity  to  the  UAV’s  air  velocity,  to  achieve  the  drag-to-lift  ratio  in  terms  of  the 
ground  velocity.  The  next  step  is  to  construct  a  generalized  result  for  costt/  and  cost-f ‘ , 

for  the  case  of  a  three-dimensional  wind  velocity.  This  is  necessary  for  our  three- 
dimensional  network  representation  of  wind,  paths  and  waypoints  in  the  various 
optimization  models. 

The  effect  of  wind  can  be  regarded  as  a  shift  in  the  frame  of  reference  of  the  UAV 
by  the  wind  velocity.  This  effect  leads  to  a  different  distance  traveled  in  the  earth’s  frame 
of  reference,  compared  to  the  one  perceived  by  the  UAV  in  its  flight  through  the  wind. 
This  study  assumes  that  flying  through  a  changing  wind  field  will  not  cause  instabilities, 
such  as  the  case  may  be  for  micro  UAVs  in  turbulence. 

The  network  will  be  divided  into  sufficiently  small  enough  arcs  to  assume  a 
constant  wind  along  an  arbitrary  arc  (i,j)  .  Let  U  be  the  UAV’s  velocity  in  Earth’s  frame 

of  reference,  Vwind  be  the  wind’s  velocity  in  Earth’s  frame  of  reference  and  URel  be  the 
UAV’s  air  velocity,  i.e.,  its  velocity  in  the  wind’s  frame  of  reference.  Therefore, 

(27)  Vm-d+VM=Vt 

The  angles/? and  y..  are  defined  as  the  direction  of  the  wind  and  the  arc  (/,/) 

respectively,  on  a  common  plane  relative  to  a  common  arbitrary  axis.  The  direction  of 
flight  in  the  wind’s  frame  of  reference  is  constrained  to  lead  to  the  next  node  j ,  and  is 
denoted  by  a .  In  other  words,  the  UAV’s  air  velocity  is  chosen  to  be  on  the  plane 
created  by  the  wind  velocity  and  the  line  that  connects  node  i  to  j ,  and  its  size  and 
direction  are  chosen  such  that  the  UAV  ends  up  at  j  .  This  is  illustrated  in  Figure  13. 
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Figure  13.  Wind,  air  and  ground  velocities  ( Vwmd ,  FRel  and  K  respectively)  and  the 
corresponding  angles  J3 ,  a  and  y,,  with  respect  to  the  horizon,  on  the  plane 
created  by  Vwind  and  VtJ .  Two  cases  are  presented  for  clarification:  (a)  Tail 
wind  (b)  Head  wind. 


The  total  travel  distance  that  the  UAV  “experiences”  in  the  frame  of  reference  of 
the  wind,  XRel  can  be  expressed  as: 

(28)  XRel=VRerty=VRel^ 


Here,  VRel  = 


V  Rel 


V^ 


Vj  are  the  airspeed  and  the  ground  speed,  respectively.  ttj  is  the 

flight  time  between  i  and  j  ,  and  Xtj  is  the  true  ground  distance  between  i  and  j  .  With 

this  in  hand,  the  travel  cost  can  be  expressed  as: 

f  D} 

^  L 


(29) 


costy  =  W  •  ri2gine  •  fRcl  ~  ^  D 

Vij 


x;; 


=  W  ■  --f- ■  [A  Vj  +  b  /  vj] 


v„ 


The  ground  speed  between  i  and  j ,  V:j ,  can  be  found  by  vector  calculations: 
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(30)  ^  =F„„  -cosOS -r„)  +  /M2 -VwJSmHp-n)  , 


where  F^  = 


F  (find 


F.  depends  only  on  the  scalar  difference  between  ft  and  y.. ,  which  justifies  the 

two  dimensional  treatment  of  the  wind  and  arc  vectors  on  the  plane  that  contains  them. 
FRel  is  the  only  controllable  variable  in  this  analysis,  and  it  will  be  chosen  to  minimize  the 
following  expression  for  the  cost: 


(31)  cost°p  =MIN< 


w -n  x  ■ V  ■ X  -\ AV  2+B/V  2 

"  '/ engine  V  Rel  ^  ij  f71  *Rel  ^  u  1  v  Rel 


Kind  •  cos(/?  -y)  +  yjvm2  -  VWin(f  sin2  (/?  -  y) 


where  the  engine’s  efficiency,  rj  ine ,  may  be  given  as  a  function  of  the  airspeed 

In  the  case  of  a  constant  flight  speed,  this  will  not  be  a  minimization  function,  and 
the  only  controlled  variable  will  be  the  direction  of  flight,  a ,  which  can  be  calculated 
after  observing  that  (see  Figure  13): 

(32)  VWind  ■  sin(/?  -  y)  =  VRel  ■  sin(y  -  a) , 


or 


(33) 


a  =  y  —  asin 


f  V 

w 

Ik, 


sin  (J3~y) 

J 


Again,  a  is  measured  relative  to  the  horizon  on  the  plane  defined  by  ft  and  ytj . 

Practically,  a  flight  direction  a  can  be  achieved  dynamically  by  the  control 
systems  of  the  UAV,  given  that  it  is  able  to  aim  to  the  next  waypoint  j  through  its 
navigation  system,  and  correct  deviations  from  the  i  -  j  line  when  needed. 


3.  Constraints  on  Airspeed 

There  are  four  constraints  on  the  UAV’s  airspeed,  which  are  developed  in  this 
section.  Any  aircraft’s  airspeed  is  limited  from  below  by  its  stall  velocity,  FStall,  under 
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which  the  lift  force  becomes  too  low  to  resist  the  aircraft’s  weight.  Hence,  our  first 
constraint  will  be:  URel  >  Fstalll .  Practically,  this  constraint  needs  to  be  further  restricted  to 

a  higher  bound  for  the  airspeed,  in  order  to  keep  the  UAV  well  on  the  safe  side,  above  its 
stall  speed. 

The  second  constraint  is  found  by  ensuring  that  the  UAV  is  able  to  overcome  the 
wind  and  fly  on  the  path  from  i  to  j  .  This  is  found  through  geometry,  as  shown  in  Figure 
14.  This  constraint  is  divided  into  two  cases: 


a) 

b) 


If  \P~Yij  | <90,  then  URel 
If  \P~Yy  | >90,  then  URel 


>K 


Wind 


sin  (P-Yil) 


>  V 

Wind  ’ 


Figure  14. 


Finding  the  lower  bounds  on  the  airspeed, 


VRel 


using  geometry,  for  two 


cases:  (a)  |  P-ytj  |  <  90  and  (b)  \f3-yiJ\>9Q). 


The  third  constraint  will  be  another  lower-bound  on  the  UAV’s  ground  speed, 
which  ensures  that  the  UAV  will  cross  any  arc  in  a  reasonable  time,  or  in  terms  of  ground 
speed:  K  >  Vmm ,  where  Vmm  is  the  lowest  allowable  ground  speed.  This  constraint  can  be 
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expressed  in  terms  of  FRel  and  Fwind  because  V ,,  is  simply  the  projection  of  V rCi  on  the 
path  i  -  j  plus  the  projection  of  V  wind  on  that  path.  Altogether: 


I'm,,  cos  >  r„„ 


or  in  terms  of  hKd : 

^Rel  —  V  ^min  _  ^Knin  Kvind  COS(/?— /)  +  Vwind 

The  fourth  constraint  arises  from  the  limitations  of  the  engine  and  takes  the  form 
FRe j  <  Vnmx ,  where  ljn,ix  is  the  maximal  achievable  airspeed  of  the  aircraft. 


The  objective  function  for  the  arc  cost  and  the  constraints  on  the  UAV’s  airspeed 
can  be  summarized  as  the  following  non-linear  programming  problem: 


c,  =  costf  =MIN 

1  1  K., 


V  >  V„  ,  >  V* 

r  max  —  r  Rel  ~  ’ 


(34) 


where : 


V*  =  < 


MAX 


MAX 


^Stalll  ’  Kvincl  sin(/^  Yij  ) 

-2VmmVWindcos (f-r)  +  vVl 


V  V 

r  Stall!  >  v  Wind  » 


2 

Wind  J 
\ 


^Vmf-2VmmVWindCos{P-y)  +  VM 


2 

Wind  J 


\P~Yij  1-90. 

P~Yj\>  90- 


4.  Turning 

The  networks  implemented  in  this  study  measure  dozens  of  kilometers,  while  the 
arcs  are  1-1.42  km  long.  It  will  be  shown  that  with  arcs  of  1  km  length  or  more,  the  cost 
of  a  turn  from  arc  (/, /)  to  arc  (j,k)  is  negligible  compared  to  the  total  cost  of  crossing 
the  arcs  (/,  j)  and  (j,k) .  To  show  that  the  cost  of  turns  is  negligible  in  the  models  of  this 
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study,  a  simplified  model  for  turns  is  considered  with  the  following  assumptions:  (a) 
turns  are  made  along  a  circular  arc,  and  (b)  turns  are  characterized  by  a  constant  radial 
acceleration  and  a  constant  absolute  speed. 


A  turn  from  arc  (/,  j)  to  arc  (j,  k)  starts  with  velocity  fj  and  ends  with 


velocity  V2 ,  such  that 


under  assumption  (b).  The  radius  of  the  turn  is  R  ,  and  its 


radial  angle#,  as  illustrated  in  Figure  15. 


j 


Figure  15.  A  circular  turn  from  arc  (i,j) to  arc  ( j,k ),  with  turn  radius  R,  over  0 
degrees. 


A  turn  is  achieved  by  producing  a  radial  acceleration  through  the  aerodynamic 
forces  created  by  the  UAV’s  wings.  In  a  two  dimensional  turn,  the  horizontal  projection 
of  the  lift  on  the  line  that  connects  the  UAV  and  the  center  of  the  desired  turn  circle 
equals  to  the  radial  force.  For  every  extra  unit  of  force  that  is  taken  from  the  lift  for  the 
purpose  of  turning,  an  extra  D/L  ■  rf^gine  units  of  thrust  need  to  be  produced  by  the  engine. 

During  a  turn,  the  UAV  must  use  a  roll  angle,  8 ,  to  create  a  projection  of  the  lift  on  the 
horizontal  plane.  Therefore,  the  extra  energy  spent  during  a  turn,  A Eijk ,  is  approximated 

as: 
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AE„k  =\hirndl  =m  DSm5  \aTumdi  =K^aTumdl  =KaTlim\dl 


(35) 


^  engine  J 

\AV\  /rAFli?©  k\AV  RQ 

-  KaTumR 0  =  k  ' - -  RQ  =  — - - =  —  - 

Twn  At  Q/co  e/(V  /R) 


=  m 


^T>sin<^ 

^  engine  J 


V  AV , 


with  the  following  definitions: 

\v\  =  \vx\  =  \v2\  and  |AF|  = 


k  =  m 


^.DsinA^ 


y  ^^1  engine  J 


Vx-V2 


The  fourth  equality  in  eq.  (35)  is  a  result  of  the  constant  acceleration  assumption,  and  the 
fifth  equality  is  a  result  of  the  circular  turn  assumption. 

The  ratio  between  A Ejjk  and  the  cost  of  crossing  arc  (/,  j)  can  be  calculated  by 
dividing  eq.  (35)  by  eq.  (25): 


(36) 


x  \MAX  m 

A E...  A 
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Vc°st  uj 


Z)sin  o 

engine  J 


mgX 
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f  D  ^ 

y  engine  J 


The  largest  angles  in  the  network  are  109.5  degrees  (e.g.,  flying  to  an  upper-right 
comer  and  returning  to  the  lower  left  corner  in  the  next  9-nodes  layer).  This  corresponds 
to  a  maximum  of  |Ak|  =  1.63|F| .  The  steepest  roll  angles  needed  for  a  turn  from  arc  to  arc 
will  be  bounded  from  above  by  8  =  30  degrees,  assuming  long  turns  (over  hundreds  of 
meters).  Therefore,  the  maximal  energy  ratio  is: 


(37) 


f  A E...  A 

Uk 

vcosW 


< 


0.82|k|: 


Since  |F|  is  usually  in  the  order  of  20-25  m/s  and  1  <  Xi}  <1.42  km,  this  ratio  will 
be  lower  than: 


0.82x25 

9.81-1000 


=  5.22% 
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in  most  cases,  and  hence,  neglected  in  this  work.  The  price  of  considering  turns  in  the 
models  of  this  study  is  compromising  the  exact  or  approximate  algorithms  in  chapter  IV, 
for  heuristic  ones.  An  exact  algorithm  that  considers  turns  would  have  required 
multiplying  nodes  in  every  location  by  the  number  of  possible  turns  to  these  nodes,  which 
may  significantly  slow  the  calculations,  and  hence,  is  not  done  in  this  work. 

5.  Climbing  and  Descending 

The  cost  of  vertical  movement  between  two  altitudes  ht  and  /?  ,  is  approximated 

as  the  original  cost  added  to  the  difference  in  the  potential  energies,  z ng[hj  -  /?,  j .  Under 

the  assumptions  that  excess  energy  cannot  be  transferred  between  arcs  due  to  the  speed 
control,  and  that  climbing  /  descending  angles  are  shallow4,  the  modified  cost  becomes: 

(38)  costfjlimb/Descend  =  costy  +  MAX { mg ( hj  -h), o}  . 

Steep  climbing  angles  may  change  the  drag  to  lift  ratio  for  a  given  speed,  which 
makes  climbing  and  descending  a  non-symmetric  processes  in  terms  of  the  extra  energy 
that  is  consumed  or  gained,  in  addition  to  the  basic  arc  cost  [29]  However,  our  use  of 
small  vertical  angles  keeps  eq.  (38)  in  a  simple  form. 

6.  Effect  of  Flight  Altitude 

In  general,  the  following  effects  can  be  identified  with  the  addition  of  the  vertical 
degree  of  freedom  for  flight: 

•  A  higher  altitude  means  a  lower  air  density,  and  a  different  drag  to  lift 
ratio. 

•  Better  wind  directions  or  speeds  may  be  found  at  higher  or  lower  altitudes, 
similarly  to  the  benefit  of  a  horizontal  degree  of  freedom.  In  general,  the 
higher  the  altitude,  the  higher  the  wind  speeds. 

The  effect  of  the  change  in  air  density  is  discussed  in  this  section. 


4  In  this  study,  the  climbing  or  descending  angles  are  bounded  by  roughly  5.7  degrees,  since  the 
vertical  layers  separation  is  100  m,  and  the  minimum  horizontal  separation  is  1000  m. 
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In  eq.  (18)  and  eq.  (19)  for  the  lumped  parameters  A  and  B ,  the  fact  that  the  air 
density  p  is  a  function  of  the  altitude  h  ,  was  stressed.  The  models  in  this  thesis  use  the 
following  approximation  in  eq.  (39)  for  the  air  density  [32]: 


(39) 


p(h)  = 


Mp(h )  _  M 
RT(h )  ~  R 


Pv 


1  + 


Lh 


lo  ) 


gM_ 

RL 


T0  +  Lh 


where: 


p0  is  the  sea  level  atmospheric  pressure  [Pa], 


T0  is  the  sea  level  standard  temperature  [K], 

L  is  the  temperature  lapse  rate  [K/m], 

R  is  the  universal  gas  constant  [J/mol/K], 

M  is  the  molecular  weight  of  dry  air  [kg/mol]. 

Use  of  the  UAV  parameters  in  Table  2,  zero  wind  and  an  arc  length  of  1  km 
yields  the  cost  function  in  Figure  16. 


Table  2.  Representative  values  for  a  small  UAV 


Parameter 

Symbol 

Value 

Mass  [kg] 

m 

5 

Wind  span  [m] 

b 

1.93 

Parasite  area  [m2] 

f 

0.028 

Oswald’s  factor 

e 

0.7 

Engine’s  efficiency 

V 

0.7 
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Arc  Cost  at  Different  Altitudes  and  Speeds 


Figure  16.  Arc  costs  at  altitudes  up  to  2000  m,  for  three  different  flight  speeds,  with  the 
UAV  parameters  given  in  Table  2. 

Figure  16  illustrates  that  the  cost  function  decreases  with  altitude  almost  linearly 
for  the  entire  relevant  range.  In  sections  IV.B  and  IV.C  this  fact  is  used  to  shorten  the 
computation  time  when  complex  algorithms  or  when  large  networks  are  involved.  This 
will  be  done  by  the  use  of  tables  and  fitting,  instead  of  repeatedly  solving  the  constrained 
nonlinear  optimization  problems  in  eq.  (34). 

E.  ESTIMATION  OF  UAV  PARAMETERS 

In  order  to  simulate  representative  UAVs  of  different  masses,  according  to  the 
different  models  (see  Chapter  IV),  the  following  approximation  for  the  wing  span  is  used: 
(40)  bn  1.041  -m°- 382  , 

where  b  is  the  wing  span  in  meters  and  m  is  the  takeoff  mass,  in  kilograms. 

Oswald’s  efficiency  factor  was  taken  to  be  0.7  and  the  parasite  surface  area,/,  is 
approximated  through  the  wing  span: 
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(41)  /□ 


f  b_  y 

ylOy  ’ 

which  approximately  coincides,  for  example,  with  Rascal  UAV  (b  □  2.4m  ,  /  0.054 nr ) 

[14]. 

Results  of  the  simulations  are  given  in  the  form  of  the  ratio  between  the  energy 
cost  of  the  suggested  paths,  compared  to  the  energy  cost  of  the  straight-line  path  with  a 
constant  flight  velocity.  Therefore  the  engine’s  efficiency  rjengine  that  appears  as  an 

inverse  factor  in  the  nonlinear  minimization  program  for  the  arc  cost  (eq.  (34))  will  have 
no  significance  for  the  results. 
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IV.  NETWORK-SCALE  OPTIMIZATION 


A.  CHAPTER  OVERVIEW 

This  chapter  discusses  the  construction  of  five  network  models  for  minimum- 
energy  paths  and  compares  their  energy  conservation  performance.  The  five  network 
models  are  (1)  a  deterministic  model  that  assumes  perfect  knowledge  of  the  atmosphere 
and  provides  an  upper  bound  on  the  optimal  energy  saving  for  comparison  with  the  other 
more  realistic  models,  (2)  a  stochastic-static  model  that  assumes  the  use  of  forecasts  only 
and  provides  a  pre-flight  plan,  (3)  a  stochastic-dynamic  model  that  takes  into  account 
real-time  measurements  by  the  GPS-INS,  (4)  an  improved  stochastic-dynamic  model  that 
considers  the  correlation  structure  of  the  atmosphere,  and  (5)  a  maximum-duration  model 
that  is  based  on  the  stochastic-dynamic  shortest  path  and  provides  an  optimized  ratio 
between  flight  duration  and  energy  consumed. 

B.  MODEL  I:  DETERMINISTIC  SHORTEST  PATH 

The  deterministic  shortest-path  model  uses  statistical  knowledge  about  the  given 
wind  forecast  to  create  a  realistic  scenario.  Dijkstra’s  shortest-path  algorithm  [31]  is  then 
used  to  find  the  path  with  the  minimum  energy  cost  from  origin  s  to  destination  t.  This 
corresponds  to  the  ideal  case,  where  complete  information  about  the  network  is  given, 
and  serves  as  an  upper  bound  on  the  optimal  energy  saving  for  the  other  models. 

The  deterministic  shortest-path  model  is  solved  after  computing  the  cost  and 
speed  to  fly  for  each  arc  in  the  network  using  the  constrained  non-linear  programming 
model  in  eq.  (34).  For  faster  calculations,  finding  the  arc  costs  may  also  be  done  through 
the  use  of  a  table  for  the  direct  translation  of  wind  velocity  to  an  arc  cost  and  a  suggested 
speed  to  fly,  based  on  the  non  linear  program  in  eq.  (34)  .  Since  the  cost  function  depends 
also  on  the  altitude,  this  study  proposes  using  a  linear  interpolation  of  two  tables,  for  two 
extreme  flight  altitudes  at  the  lower  and  upper  potential  flight  levels.  The  example  in 
Figure  16  demonstrates  that  a  linear  behavior  will  be  a  good  approximation  between  sea 
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level  and  2000  m.  It  should  be  noted  that  different  tables  must  be  constructed  for  each 
UAV,  as  the  cost  function  depends  on  its  physical  parameters.  An  example  for  such 
tables  is  shown  in  Figure  17. 


Arc  Cost  J 


5  10  15  20  25 

Mean  Wind  Speed  [m/s] 


Speed  to  Fly  m/s 


5  10  15  20  25 

Mean  Wind  Speed  [m/s] 


(a) 


(b) 


Figure  17.  (a)  Arc  costs  [J]  and  (b)  Suggested  speed  to  fly  as  functions  of  wind  speed 

and  relative  wind  direction  for  the  UAV  parameters  in  Table  1. 


The  deterministic  shortest-path  model  uses  eq.  (10)  for  one  realization  of  the 
distribution  of  wind  velocities,  given  the  forecasted  components  of  the  wind  velocity.  For 
consistency  when  comparing  this  model  with  the  other  models,  the  same  realizations  of 
wind  velocities  are  used  in  order  to  provide  a  valid  upper  bound  by  model  I. 

The  steps  for  implementing  model  I  are  as  follows: 

1.  Construct  a  network  for  the  flight  space,  with  origin  node  seN  and 
destination  node  teN,  according  to  section  III.B. 

2.  Interpolate  a  forecast  in  the  spatial  and  temporal  domain  limits  according 
to  section  II.E. 

3.  Create  a  realization  of  wind  velocities,  using  eq.  (10). 

4.  For  each  arc  in  the  network,  calculate  its  cost,  according  to  eq.  (34). 

5.  Calculate  the  shortest  path  from  s  to  t  using  Dijkstra’s  algorithm  [31]. 

C.  MODEL  II:  STOCHASTIC-STATIC  SHORTEST  PATH 

The  stochastic-static  shortest-path  model  uses  the  statistical  behavior  of  the  wind, 
to  calculate  the  expected  (average)  arc  costs  in  the  network,  and  use  these  to  find  the 
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shortest  path  with  Dijkstra’s  shortest  path  algorithm.  This  model  represents  the  case  of  a 
pre-flight  plan,  and  assumes  no  updates  during  the  flight.  It  serves  as  a  lower  bound  on 
the  optimal  energy  saving  for  the  stochastic-dynamic  model,  as  it  uses  only  forecast  data 
and  no  updates  of  the  path  during  the  flight. 

The  expected  arc  cost  is  expressed  through  the  following  equation: 


(42) 


E[costv]  = 

00  ^  ^  V 

=  f  f  costij[vwind,(p\Vwind(p  +  p0),xAg{p)f(VWind)dpdVWind, 


VtVind  P—7U 


where  cost.,  is  the  arc  cost,  as  in  eq.  (34),  given  as  a  function  of  the  wind  speed 


VWmd  =  |  V  wind  I  and  its  relative  direction  to  the  arc  (i,j)  vector  X ,, .  j80  is  the  average 
direction  of  the  wind,  and  ft  is  the  fluctuation  around  it,  that  follows  a  normal 
distribution  gift)  as  in  eq.  (9).  <p  is  the  relative  angle  between  Xy  and  the  fluctuating 
wind  velocity,  V wm  .  The  wind  speed  follows  a  Weibull  distribution  f(VWind),  as  in  eq. 
(1). 


For  strong  enough  head  winds,  no  feasible  solution  could  be  found  to  the 
constrained  nonlinear  program,  since  a  minimum  ground  speed  would  not  be  achieved. 
Therefore,  the  integration  limits  are  reduced  into  the  range  of  two  standard  deviations 
below  and  above  the  average  wind  speed  and  direction,  using  av  and  ap  respectively,  as 
shown  in  eq.  (43).  For  the  wind  direction  [J ,  the  distribution  is  nonnal  and  hence,  this 
approximation  “skips”  4.55%  of  the  range  of  angles  in  the  integration,  for  a  range  of  two 
standard  deviations  around  the  mean  (i.e.,  there  is  a  4.55%  chance  to  encounter  wind 
directions  outside  the  range  of  the  integral).  The  wind  speed  follows  Weibull  distribution 
which  has  the  cumulative  distribution  in  eq.  (2).  Figure  18  shows  that  the  resulting 
reduction  in  the  range  of  integration  of  the  approximation  is  quite  low,  in  the  order  of 
4.5%.  The  relation  between  the  average  wind  speed  and  its  standard  deviation  is  found 
through  eq.  (10). 
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E\costq\  = 

(43)  V Wind +2 (Jv  Po  +2(7 P 

=  J  J  cosUj^wind^^d{P^ P^,Xi^g{p)f{ymnd)dpdVWind. 


V Wind  —2(Jy  J3q  2<J p 


where  V wm  is  the  average  wind  speed. 


Approximation  Relative  to  Full  Integration  on  Windspeed 


Figure  18.  Comparison  of  approximation  (43)  to  the  full  integral  in  eq.  (42),  regarding 
the  range  of  integration  over  wind  speeds. 

The  expressions  within  the  integrals  may  produce  infeasible  solutions.  In  order  to 
avoid  averaging  of  solutions  that  include  infeasible  airspeeds,  the  constraints  of  eq.  (34) 
are  first  ignored.  Then,  whenever  the  solution  for  the  unconstrained  nonlinear  program  is 
higher  than  the  maximum  allowed  airspeed,  the  maximum  airspeeds  will  be  used  as  the 
speed  to  fly,  and  the  high  cost  functions  will  serve  as  penalties.  That  way,  arcs  that  can 
usually  be  crossed  with  a  legitimate  flight  speed  (under  non-extreme  wind  velocities)  will 
remain  available  on  one  hand,  but  penalized  with  high  cost  on  the  other  hand.  Figure  19 
illustrates  such  a  situation,  where  flight  speeds  are  bounded  by  38  m/s,  and  the 
corresponding  cost  of  solutions  that  may  require  speeds  higher  than  that  are  penalized. 
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Since  eq.  (34)  solves  a  non-linear  program  for  each  value  of  /?  and  VWind  in  the 

range  of  the  double  integral,  it  will  be  computationally  expensive  to  solve  it  in  real-time 
for  every  arc  in  the  network.  There  are  two  efficient  ways  of  calculating  the  expected 
cost: 

1.  Draw  enough  random  realizations  from  g(j3 )  and  /  (VWind  )  to  estimate  the 
average. 

2.  Construct  a  table  for  the  average  cost  for  every  combination  of  VWind  and 
cp ,  and  use  a  linear  interpolation  for  the  dependency  on  the  flight  altitude, 
as  described  in  Model  I. 

The  first  method  will  be  used  when  comparing  the  three  models,  due  to  the  fact 
that  a  set  of  realizations  for  the  random  fluctuations  of  the  wind  (defined  by  VWjnd  and  (p  ) 
are  already  available  from  model  III  (see  next  section).  This  study  suggests  using  the 
second  method  in  real-time  calculations. 

The  suggested  speed  to  fly  will  be  the  expected  best  speed  to  fly,  for  the  entire 
range  of  realizations  from  gift)  and  /  (  VWind  ) .  This  may  not  be  the  optimal  choice  for  the 

flight  speed,  but  serves  as  a  representative  value  for  it.  In  this  model,  no  information  on 
the  real  wind  velocity  will  be  measured  during  the  flight.  Thus,  the  exact  optimal  flight 
speed  is  unknown.  Figure  19  shows  the  results  for  the  average  arc  cost  and  suggested 
speeds  to  fly  for  a  specific  UAV. 


Figure  19.  (a)  Expected  arc  costs  [J]  and  (b)  Expected  speed  to  fly  as  functions  of  wind 

speed  and  relative  wind  direction  for  the  UAV  parameters  in  Table  1. 
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The  steps  of  implementing  model  II  are  as  follows: 

1.  Construct  a  network  for  the  flight  space,  with  origin  S  and  destination  T, 

according  to  section  III.B. 

2.  Interpolate  a  forecast  according  to  section  II.E., 

3.  Option  A: 

1)  Create  a  large  set  of  realization  of  wind  velocities,  using  eq.  (10). 

2)  For  each  arc  in  the  network  calculate  its  set  of  realized  costs  according 

to  eq.  (34)  and  take  their  average  to  obtain  the  expected  arc  costs. 

3)  Penalize  costs  for  unfeasible  solutions,  as  shown  in  Figure  19. 

Option  B: 

1)  Create  a  table  of  expected  arc  costs  and  speeds  to  fly,  using  eq.  (43). 

2)  Penalize  costs  for  unfeasible  solutions,  as  shown  in  Figure  19. 

4.  Calculate  the  shortest  path  from  s  to  t  using  the  expected  arc  costs  and 

Dijkstra’s  algorithm  [31], 

D.  MODEL  III:  STOCHASTIC-DYNAMIC  SHORTEST  PATH 

The  stochastic-dynamic  shortest-path  model  uses  an  iterative  approximate 
dynamic  programming  algorithm  that  provides  the  expected  path  cost  from  each  node  to 
the  destination  node.  This  algorithm  dynamically  re-optimizes  the  UAV’s  route  based  on 
dynamically  collected  information  about  the  local  wind  field.  The  algorithm  used  is  based 
on  [32]  (See  also  [34]).  Arc  costs  are  given  as  random  variables  based  on  the  distributions 
of  the  wind  velocities.  When  UAV  reaches  a  node,  it  measures  the  local  wind  realization, 
and  infers  the  cost  of  each  of  the  neighboring  arcs  (in  its  forward  star).  Then,  it  makes  a 
decision  where  to  go  next  based  on  the  measured  arc  costs  and  the  expected  cost-to-go 
(cost  until  the  destination  node)  for  each  of  the  neighboring  nodes.  The  expected  cost-to- 
go  on  node  i  is  denoted  by  V .  Bellman’s  equation  [33]  gives  the  optimality  condition 
for  the  V 's  .  In  terms  of  our  network  formulation,  Bellman’s  equation  translates  into: 

(44)  Vj=E  min  (vJ '  +  cost \.)  , 
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where  j  e  {/ 1  (i,l)  e  A]  are  the  neighbors  of  node  i  in  its  forward  star  [33].  Here,  the 

expectation  E  is  used  since  the  arc  cost  is  a  random  variable  that  has  an  infinite  number 
of  potential  realizations.  Therefore,  the  model  provides  the  path  with  the  minimum 
expected  energy  cost. 

Solving  V  for  every  node  in  the  network  would  be  possible  in  a  deterministic 

scenario,  using  recursion  from  the  destination  node  backwards.  In  the  case  of  discrete 
distribution  for  the  arc  costs,  it  may  also  be  possible  to  solve  eq.  (44)  exactly  by 
recursion.  In  our  case,  the  cost  function  has  a  continuous  distribution,  and  an  exact 
solution  for  eq.  (44)  cannot  be  computed  in  finite  computing  time.  Therefore,  this  study 
uses  a  Monte-Carlo  based  algorithm  with  the  converging  cost-to-to  estimates  in  eq.  (46), 
see  below,  to  replace  eq.  (44).  The  expected  cost-to-go  is  iteratively  updated  with  a 
randomly  generated  realization  of  the  on-board  wind  velocity  measurement.  Iterations 
terminate  upon  reaching  a  reasonable  approximation  of  the  expected  cost-to-go.  Two 
examples  for  the  convergence  of  V ,  where  i  is  the  origin  node,  are  plotted  in  Figure  20. 
In  both  cases,  the  initial  error  is  no  more  than  0.5%  of  the  total  cost.  Convergence  into  the 
range  of  0.03%  is  achieved  after  20  iterations. 


Figure  20.  The  expected  cost-to-go  convergence  of  the  origin  nodes  in  two  sets  of 
network  realizations. 
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Measurement  errors  in  the  dynamic  estimation  of  the  local  wind  velocity  are 
simulated  through  the  use  of  available  data  in  the  literature  regarding  the  accuracy  of 
GPS-INS  systems  [37]-[40]  and  estimations  for  the  case  where  compass  corrections  for 
the  UAV’s  “crab  angle”  are  available  [14]. 

Another  type  of  error  that  is  taken  into  account  is  the  “Spatial  fluctuation  error.”  It 
represents  the  different  wind  velocities  at  different  points  along  the  arcs  emanating  from 
the  node  where  the  measurement  was  made.  The  error  is  calculated  as  the  variation  in 
wind  velocities  for  the  collection  of  all  adjacent  node  pairs  in  the  forecast,  including  on 
the  diagonals.  The  histograms  of  the  variations  in  speed  and  direction  between 
neighboring  nodes  are  plotted  in  Figure  21.  These  variations  are  added  as  errors  to  the 
measurements  by  the  GPS-INS  system  as  they  represent  the  difference  between  local 
wind  velocities  and  the  wind  velocities  along  neighboring  arcs. 


Wind  Speed  Difference  [m/s]  Wind  Direction  Difference  [degrees] 


(a) 


(b) 


Figure  21.  Wind  variation  between  adjacent  nodes  in  terms  of  (a)  wind  speed  and  (b) 
wind  direction. 


Errors  in  the  estimation  of  wind  speed  and  direction  are  assumed  to  have  a 
Normal  distribution  with  a  standard  deviation  that  equals  to  the  root  sum  squares  (RSS) 
of  the  two  types: 
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(45) 


cr 


=  a  cr, 


GPS-INS 


+  cr 


Fluctuations 


where  i  stands  for  either  speed  or  direction,  <Jestimate  is  the  total  estimation  error, 

& gps-ins  >s  ^e  standard  deviation  of  the  GPS-INS  error,  and  o fluctuations  is  the  standard 

deviation  of  the  spatial  fluctuations.  Table  3  summarizes  the  data  on  measurement  errors 
and  spatial  fluctuations,  collected  from  the  Marina  weather  station  and  COAMPS 
forecasts. 


Table  3.  Estimation  errors  of  wind  speed  and  direction  in  the  fonn  of  standard 
deviations 


Parameter  (standard  deviations) 

Symbol 

Value 

Speed  Measurement  Error  [m/s] 

speed 

G GPS-INS 

2 

Speed  Spatial  Fluctuation  Error  [m/s] 

speed 

fluctuation 

0.7 

Total  Speed  Error  [m/s] 

speed 
rj  . 

estimate 

2.1 

Direction  Measurement  Error  [degrees] 

direction 

G GPS-INS 

10 

Direction  Spatial  Fluctuation  Error  [degrees] 

direction 

fluctuation 

5 

Total  Direction  Error  [degrees] 

direction 
^  estimate 

11.2 

The  steps  of  implementing  model  III  are  as  follows: 

1.  Construct  a  network  for  the  flight  space,  with  origin  s  and  destination  t, 
according  to  section  III.B. 

2.  Interpolate  a  forecast  according  to  section  II.E. 

3.  Start  with  an  initial  guess  for  the  expected  cost-to-end,  V' ,  for  every  node 
i  in  the  network,  e.g.,  using  the  results  of  model  II  to  find  the  cost  of  the 
shortest  path  over  the  expected  values  of  the  arc  costs,  from  each  node  i 
to  the  destination  t. 

4.  Create  a  large  set  of  n  »  1  realization  of  wind  velocities,  using  eq.  (10). 

5.  For  each  realization,  with  index  ke{l, calculate  the  cost  of  every 
arc  (i,j)eA,  cost* ,  according  to  eq.  (34)  using  a  sample  of  the  wind 
velocity  on  arc  (/,  j)  . 

6.  For  every  node  i  in  the  network,  update  the  value  of  V,'  ,  to  Vlk  as 
follows: 
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(46)  V'k={\~  ak )  Min\  If.  +  cost*  ]  +  akVf , 

where  ak  serves  as  a  weighting  factor,  that  can  be  computed  in  various  ways.  To  achieve 
an  arithmetic  mean  at  every  iteration  k ,  use: 


(47) 

7. 


1 

ak  = - . 

k  + 1 

Repeat  stage  6  for  the  n  realizations.  Figure  20  illustrates  the  rate  of 
convergence  of  Vs . 


Now  the  model  is  ready  for  a  flight  with  dynamic  measurements: 

1 .  Start  at  node  5. 

2.  Measure  the  real  arcs  costs  cost™1  from  the  current  node  i  to  the  adjacent 
nodes  j  ,  for (i,  j)  e  A. 

3.  Fly  to  node  /  that  satisfies: 

(48)  /  =  Argmin  [vj  +  cost™1  J  . 

j 

4.  Repeat  stages  2-3  until  reaching  the  destination  node  t. 

In  a  simulation,  stages  2-3  use  arc  costs  with  random  measurement  errors  according  to 
Table  3. 


E.  MODEL  IV:  CORRELATION-BASED  STOCHASTIC-DYNAMIC 

SHORTEST  PATH 

This  model  uses  the  spatial  and  temporal  correlation  structure  of  the  atmosphere, 
as  described  in  section  II.F,  to  achieve  reduced  energy  consumption.  This  model  assumes 
a  spatially-autocorrelated  joint  distribution  for  the  wind  velocities,  instead  of  independent 
marginal  distributions  as  was  assumed  in  the  previous  models.  This  model  relates  only  to 
the  correlation  structure  of  the  fluctuations  around  the  forecasted  quantities.  The  model  is 
similar  to  model  IV,  with  the  difference  that  measurements  give  us  new  information 
about  wind  in  the  entire  space  of  flight,  and  not  only  on  the  adjacent  arcs.  As  noted  in 
section  II.F,  the  temporal  and  spatial  correlation  structure  of  the  atmosphere  is  not 
available  at  the  time  of  the  study,  not  even  for  a  specific  region,  and  may  require 
experimentation  or  construction  of  a  physical  theory. 


44 


The  model  uses  techniques  for  creating  correlated  wind  velocities,  through 
correlated  normal  and  Weibull  distributions  (see  [41], [42]  and  Appendix)  for  the  wind 
direction  and  speed,  given  a  correlation  structure  of  the  wind  for  the  relevant  ranges  of 
space  and  time  (see  section  ILF  for  examples  of  one-dimensional  correlation  structures). 

The  steps  of  implementing  model  IV  are  as  follows: 

1.  Construct  a  network  for  the  flight  space,  with  origin  s  and  destination  t, 
according  to  section  III.B. 

2.  Interpolate  a  forecast  according  to  section  II.E., 

3.  Start  with  an  initial  guess  for  the  nodes-costs-to-end,  V0' ,  for  every  node  i 
in  the  network. 

4.  Start  at  node  s. 

5.  Create  a  large  set  of  n  realization  of  correlated  wind  velocities  from 
current  node  onwards,  towards  the  destination,  using  eq.  (10).  Create 
correlated  wind  speeds  and  correlated  wind  directions  separately, 
according  to  the  Appendix.  Then  combine  them  into  a  general  wind 
velocity. 

6.  For  each  realization,  with  index  k  g  {l,...,n} ,  calculate  the  cost  of  every 
arc  (/,  j)  e  A  ,  cost * ,  according  to  (34). 

7.  Choose  the  largest  subset  of  the  realizations  p  g  P  a  {l...n}  such  that  the 

costs  of  the  neighboring  arcs  to  the  current  node  are  within  a  pre-defined 
environment  s  of  the  corresponding  measured  arcs  costs,  e.g.,  within 
£■=5%  5. 

8.  Using  the  subset  of  realizations  in  7,  update  the  value  of  V'  according  to 
(46)  for  every  node  i  between  the  current  node  and  the  destination,  over 
|P|  iterations,  as  described  in  model  III. 

9.  After  updating  the  V’  ’s,  choose  the  next  node  to  fly  to,  according  to  (48) 

and  the  measurements,  and  repeat  stages  7  and  8  until  reaching  the 
destination  t. 


5  In  a  high  resolution  network,  all  of  the  nine  neighboring  arcs  should  have  a  similar  wind  velocity, 
both  in  the  measurements  and  in  the  created  set  of  realizations,  due  to  the  high  correlation  of  wind  velocity 
between  adjacent  arcs. 
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F. 


MODEL  V:  MAXIMUM  DURATION  STOCHASTIC-DYNAMIC 
SHORTEST  PATH 


The  maximum  duration  stochastic-dynamic  shortest  path  model  is  based  on  the 
stochastic-dynamic  model.  It  also  assumes  a  flight  between  origin  s  and  destination  t. 
Instead  of  minimizing  energy  per  distance,  it  maximizes  the  flight  time  per  unit  energy 
(e.g.,  hours  of  flight  per  gallon  of  fuel,  or  minimizes  energy  per  unit  time).  This  model 
follows  the  steps  of  model  III  exactly,  except  for  one  difference:  the  energy  cost  function 
in  eq.  (34)  is  divided  by  4f. /If  ,  to  obtain  energy  per  unit  time: 

cost?  =MIN{w-n-^{AVj+BIV^,, 

^Rel  V  ' 

(49)  subject  to: 


where  V  *  is  defined  as  in  eq.  (34). 

Let  Ej  be  the  energy  cost  of  crossing  arc  (/,/)  and  71.  be  the  time  to  cross  arc 
(/,/) .  Define  A  and  A  be  as: 


(50) 


A  =  Min 

all  paths 


z  a 

(i,j)epath 


Z  r, 

(i,j)Gpath 


and: 


(51)  A  =  Min  Y 

all  paths  , .  ,  T 

(i, j)s path  ±ij 

Solving  eq.  (49)  and  applying  a  shortest-path  algorithm  will  not  yield  an  exact  solution, 
or  a  solution  that  converges  to  the  exact  one,  as  in  models  I-IV.  Instead  of  providing  the 
path  with  the  minimum  energy  per  unit  time,  A ,  by  solving  eq.  (50),  the  algorithm  solves 
the  approximate  expression  in  eq.  (51),  to  find  the  path  with  the  minimum  sum  of  the 

ratios  E..  / 71.,  A.  Hence,  this  is  a  heuristic  algorithm  that  doesn’t  necessarily  provide  the 
overall  best  achievable  ratio  of  energy  per  unit  time.  This  arises  from  the  fact  that  the 
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energy  to  reach  from  origin  to  destination  can  be  directly  translated  into  energy  per  unit 
distance,  but  not  into  energy  per  unit  time,  due  to  the  space-based  architecture  of  the 
network.  Despite  that,  this  is  a  heuristic  algorithm  it  still  produces  a  significant  increase 
in  the  flight  duration  per  fixed  energy. 

G.  RESULTS 

The  first  three  models  were  implemented  in  MATLAB,  using  COAMPS  weather 
forecasts  over  a  160  km  by  160  km  region  at  Yucca  Air  Force  Test  Site,  Nevada.  A  total 
of  1000  North-South  and  East-West  paths  of  40  km  round  trips  (approximately  20  km 
ingress  and  20  km  egress)  were  simulated  through  a  Monte  Carlo  numerical  simulation. 
Horizontal  movement  was  constrained  by  the  network  architecture  to  5  km,  and  vertical 
movement  was  constrained  to  be  within  500  m.  Origins  and  destinations  were  chosen 
randomly  under  the  constraint  that  they  laid  no  less  than  100  m  above  the  highest  terrain 
features  in  the  area  of  flight.  The  resulting  network  contains  1120  nodes  and  5980  arcs. 
Each  of  the  models  produced  a  suggested  path  to  fly  and  an  associated  energy  cost,  which 
was  compared  to  the  energy  cost  of  a  straight-line  flight. 

In  the  following  discussion,  all  energy  savings  are  computed  relative  to  energy 
consumption  of  a  straight  line  flight  path  from  origin  to  destination  with  a  constant  flight 
speed.  Significant  energy  savings  were  achieved  in  head-winds  scenarios,  in  which  the 
models  produced  complex  paths  that  tried  to  avoid  unfavorable  strong  currents. 
Conversely,  the  models  did  not  produce  significant  energy  savings  in  tail-wind  scenarios, 
as  the  suggested  paths  were  similar  to  the  straight-line  paths.  An  example  for  these 
observations  is  shown  in  Figure  22.  The  forecasted  horizontal  component  of  the  wind 
velocity,  Vx ,  is  shown  as  a  contour  plot  in  the  background. 

Figure  23  compares  the  energy  savings  achieved  by  the  deterministic,  stochastic- 
static  and  stochastic-dynamic  models,  as  a  function  of  the  average  wind  speed  along  the 
direct  s-t  path.  As  expected,  the  deterministic  model  provides  an  upper  bound  on  the 
energy  saving  for  the  other  models.  The  stochastic-dynamic  model  produces  energy 
savings  as  good  as  the  deterministic  model,  within  less  than  10%  tolerance.  This 

performance  is  significantly  better  than  the  stochastic-static  model.  The  latter  model 
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produced  energy  savings  of  no  more  than  10%.  For  a  100  kg  UAV,  the  stochastic- 
dynamic  model  produced  an  average  improvement  (i.e.,  reduction)  of  15.1%  in  the 
energy  consumption  compared  to  the  straight-line  paths,  and  an  average  improvement  of 
30.1%  in  windy  cases,  in  which  the  average  wind  speeds  were  higher  than  15  m/s.  In 
general,  the  higher  the  wind  speeds  the  better  the  energy  conservation  that  the  models 
produce,  compared  to  the  straight-line  flight. 
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Figure  22.  Illustration  of  a  100  km  long  round  trip  from  top  and  side  views  of  the 
suggested  path  to  fly,  using  model  III. 

Figure  24  shows  the  raw  results  for  the  stochastic-dynamic  model  (model  III) 
without  averaging  them  for  similar  wind  speeds  as  in  Figure  23.  It  is  observed  that  the 
variation  in  the  results  around  the  average  values  is  significantly  higher  for  a  10  kg  UAV 
compared  to  a  100  kg  UAV.  This  translates  to  a  higher  predictability  in  the  case  of 
heavier  UAVs. 

The  maximum  duration  model  was  simulated  without  the  horizontal  degrees  of 
freedom,  to  represent  a  border-patrol-like  scenario,  in  which  the  flight  must  follow  a  long 
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and  narrow  area.  Results  for  1000  flights  are  shown  in  Figure  25.  Best  results  are 
achieved  for  wind  speeds  between  10  m/s  and  15  m/s,  in  the  range  of  20-80%  increase  in 
the  flight  duration,  compared  to  a  non-optimized  flight. 


Average  Wind  Speed  on  s-t  Path  [m/s] 


(a) 


Average  Wind  Speed  ons-f  Path  [m/s] 


Figure  23.  Averaged  results  from  1000  flight  simulations,  comparing  the  energy 
consumption  of  paths  suggested  by  models  I,  II  and  III  to  the  straight-line 
path  with  constant  speed,  as  a  function  of  the  average  wind  speed  along  the 
direct  s-t  path,  for  two  UAV  masses:  (a)  10  kg  and  (b)  100  kg. 
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Figure  24.  Results  from  1000  flight  simulations,  comparing  the  energy  consumption  of 
model  III  (Stochastic-Dynamic  shortest  path)  to  the  straight  line-path  with 
constant  speed,  as  a  function  of  the  average  wind  speed  along  the  s-t  path, 
for  two  UAV  masses:  (a)  10  kg  and  (b)  100  kg. 
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Figure  25.  Results  from  1000  flight  simulations,  comparing  the  increased  flight 
duration  per  unit  energy  of  model  V  (Stochastic-Dynamic  Maximum 
Duration  path)  to  the  straight  line-path  with  constant  speed,  as  a  function  of 
the  average  wind  speed  along  the  s-t  path,  for  a  100  kg  UAV. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

The  goal  of  this  thesis  is  to  provide  a  quantitative  analysis  of  the  potential  energy 
savings  for  unmanned  aerial  vehicles  (UAVs)  through  the  use  of  wind  currents.  This 
thesis  has  demonstrated  a  significant  average  improvement  of  30%  in  energy 
consumption  and  flight  endurance  of  UAVs,  through  network  optimization  of  flight  paths 
and  the  use  of  mesoscale  weather  forecasts  with  dynamic  measurements. 

Five  optimization  models  were  constructed-four  for  the  minimum-energy  path 
and  one  for  the  maximum  duration  path.  The  minimum-energy  models  produced 
significant  energy  savings,  especially  when  strong  head  winds  exist  and  should  be 
avoided.  When  strong  tail  winds  exist,  the  optimal  flight  path  is  usually  very  similar  to 
the  straight-line  path  and  the  models  do  not  produce  any  significant  improvement.  Pre¬ 
flight  path  plans,  based  on  weather  forecasts,  were  shown  to  be  inadequate  for  achieving 
significant  energy  saving.  With  the  use  of  dynamic  measurements  of  wind  velocities 
throughout  the  flight,  the  model  produces  significantly  better  results  that  are  much  closer 
to  the  upper  bound  of  energy  saving  given  complete  wind  knowledge.  Therefore,  a 
stochastic-dynamic  shortest-path  model  is  recommended  for  implementation  in  UAV 
flights,  both  in  operations  and  training.  This  model  will  realize  the  triple  benefits  of  (1) 
conservation  of  fuel  and  battery,  (2)  achievement  of  longer  flight  distance  and  duration, 
and  (3)  improved  operational  flexibility. 

B.  SUGGESTED  WORK  AHEAD 

There  are  several  interesting  scopes  to  be  extended  from  this  thesis.  The  first 
immediate  step  would  be  to  implement  and  verify  the  models  and  algorithms  on  an  actual 
UAV  platform.  This  can  be  achieved  through  collaboration  between  NPS  and  the  U.S. 
Naval  Research  Laboratory.  This  work  item  will  require  modification  of  existing  control 
software  in  order  to  implement  Model  III:  Stochastic  Dynamic  shortest-path. 
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The  second  work  item  would  be  to  expand  the  existing  models  to  handle  more 
complex  flight  patterns,  in  addition  to  the  simple  origin-destination  operations  assumed  in 
this  thesis. 

The  third  work  item  would  be  to  collect  data  on  the  spatial  and  temporal 
correlation  structure  of  the  wind.  Such  data  collection  is  required  for  the  implementation 
of  Model  IV:  Correlation  based  Stochastic  Dynamic  shortest-path.  This  model  is 
expected  to  provide  better  performance  than  the  Model  III:  Stochastic  Dynamic  shortest- 
path. 

The  fourth  work  item  would  be  to  conduct  sensitivity  analysis  of  the  algorithm 
performance  with  respect  to  aerodynamics  and  environmental  conditions.  Aerodynamic 
variables  should  include  mass,  wing  span,  parasite  surface  area  and  Oswald’s  efficiency 
factor.  Environmental  conditions  should  include  flight  distance,  flight  altitude,  terrain 
features  (hilly  versus  flat),  network  size  and  discretization  level.  The  last  two  conditions 
contribute  to  the  vertical  and  horizontal  degrees  of  freedom  for  flight.  However,  they 
increase  the  computation  time  of  the  algorithms  and  hence  should  be  balanced 
appropriately. 
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APPENDIX. 


Creation  of  approximately  correlated  random  variables  (RVs)  with  Normal  and 
Weibull  distributions  given  their  correlation  structure,  can  be  done  in  the  following  ways 
(based  on  [41]  for  Normal  RVs  and  on  [42]  for  Weibull  RVs). 

Creation  of  Correlated  RVs  with  Normal  Distribution: 

Let  X  =  (Xj,...,XH)  □  MFV(0,E)  be  a  random  vector  with  a  multi-variate  normal 
distribution  (MVN)  with  mean  0  and  covariance  matrix  I .  Since  X  is  symmetric  positive 
definite,  it  can  be  expressed  using  Cholesky  decomposition  as:  E  =  CTC ,  where 
C  =  \[dU  such  that  D  is  a  diagonal  matrix  with  positive  elements,  and  U  is  an  upper 
triangular  matrix.  The  way  to  generate  X  =  (V,,...,  Xi:)  □  MVN(M ,E)  is  as  follows: 

1.  Generate  a  vector  Y  =  (Yx,...,Yn)  □  MLV(0,7)  of  n  independent  RVs,  with  mean  0 
and  standard  deviation  1 . 

2.  Find  C ,  the  Cholesky  decomposition  matrix  of  E  . 

3.  X  =  M  +  CtY  is  the  required  set  of  correlated  RVs. 

Creation  of  Correlated  RVs  with  Weibull  Distribution  (  Approximation): 

The  way  to  generate  X  =  (Xl,...,Xn)  □  MVW{p,a,p) ,  where  p  =  (jul,...,jun) are  the 
averages,  cr  =  (cr,,...,cr,I)are  the  standard  deviations  and  pis  the  n  x  n  correlation  matrix 
is  as  follows: 

1 .  Define  an  n  x  1  vector  V  ,  such  that:  Vt  =  — 

Mi 

2.  Define  the  matrix  p*  through  its  elements: 

p]  =  1 .063  -  0.004py  - O.OOlp. 2  - 0.2 (jf  +  V} )  +  0.337 (Vr  +  Vf )  + 

+ 0.007  plJ(vi+ Vj )  -  o.oovvyj 
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3.  Define  ^  =  p* p 

4.  Create  Y  =  (f,,...,  f)  □  MIW(0,^»)as  described  above  for  correlated  Nonnal  RVs. 

-1.086 

,  c,  = — -r^1 — - ,  as  defined  for  the  two  parameters  of  Weibull 

r  i+— 

l  ki) 

distribution  in  equations  (6)  and  (7). 

6.  W  ={-log[l-©(i;)]}Vc,  is  the  required  multi  variate  vector  of  correlated  RVs, 
where  the  cumulative  Nonnal  distribution  ©(^defined  as: 


(  \ 

5.  Let:  kt  =  — 

UJ 
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